本推文来源：城市学派
提醒：该推文中关于中国地图存在问题，大家要学会甄别。规范表达中国地图。
引用格式： Nelson, A., Weiss, D. J., van Etten, J., Cattaneo, A., McMenomy, T. S., & Koo, J. (2019). A suite of global accessibility indicators. Scientific data, 6(1), 266.
摘要
资源和机会的充足获取对于可持续发展至关重要。提高获取，尤其是在农村地区，需要采取有效措施来改善当前对这些资源和机会所在位置的访问。近期的工作已经开发了一个全球城市地图，这些城市在 2015 年的人口超过 5 万人。然而，资源和机会的提供将在从小镇到超大城市等广泛范围的定居点中有所不同，对这些不同规模定居点的访问也应进行测量。在这里，我们呈现了 2015 年的九个全球旅行时间可达性指标，以大约 1 公里的空间分辨率，适用于一系列定居点大小类别。我们通过与 Google 驾驶方向应用中的旅程时间进行了验证，覆盖了 1,511 个 2°×2°的瓷砖，代表了 47,812 个旅程。我们观察到非常好的一致性，尽管我们的估计通常比 Google 应用的要短，中位数差异为-13.7 分钟，中位数百分比差异为-16.9%。
背景与概述
资源、服务和机会的获取，尤其是在城市中的集中，是衡量农村发展、农业生产力、食品消费者和食品生产者进入市场的机会以及贸易的重要且常用的指标 1 、 2 、 3,4 和 5 。更广泛地说，城市以及连接它们的交通网络是必不可少的基础设施，提供了人们和产品从 A 点到 B 点旅行的手段，从而促进了社会和经济互动，并提供了诸如教育和医疗保健等基本服务的交付。通常，拥有良好接入的群体具有更大的社会和经济发展的机会，较低的成本和更高的互动水平，而较差的接入则意味着更高的成本、较少的机会和较差的健康和教育结果。
不平等的获取机会可能导致更大的社会和经济差距 6 。运输网络的不恰当扩展也可能破坏自然环境，导致森林砍伐和过度开采易于获取的自然资源 7 。另一方面，合理规划的获取改善可以带来更好的农村健康、财富和经济生计结果，同时限制对环境的影响 8 。
2018 年，韦斯等人发表了关于 2015 年（联合国可持续发展目标的基准年）全球到达 5 万人或以上人口城市的时间数据集，并展示了农村福祉与到达这些城市的时间之间的负相关关系。这是对 2000 年（千年发展目标的基准年）数据集的更新和改进。2015 年的数据集利用了以下改进：描述人类聚居地大小和位置的改进数据层；更全面的交通网络和旅行速度信息；描述非道路速度的改进环境层，以及计算旅行时间时能够考虑等角矩形（经度-纬度）网格上的距离失真所需的计算工具。
2000 年和 2015 年的数据集均使用旅行时间作为易于理解的度量标准，来表示对人类聚居地的物理可达性。然而，重要的是要区分物理、经济和社交的可达性，并认识到良好的物理可达性并不意味着经济和社会可达性也同样良好。在经济可达性的情况下，一个有良好物理可达性到人类聚居地的人可能没有足够的经济手段来有效地使用交通网络。同样，同一户家庭的成员在聚居地可能面临不同的资源和机会的社交可达性水平，如果他们因文化规范而被阻止或阻止使用这些资源和机会。在这里，我们使用相同的旅行时间指标 9,10 作为物理可达性的衡量标准。
这里描述的数据在两个方面进行了进一步的改进。首先，它考虑了城市规模在资源和服务提供方面的重要性。2000 年和 2015 年数据集的一个限制是，为所有 5 万人或以上人口的城市估计了到最近城市的通勤时间。最近的城市可能是人口超过 1000 万的超大城市，也可能是 5 万人左右的小区域城市，或者两者之间的任何城市。在这个城市规模的范围内，资源、服务和机会的可用性存在巨大的差异。其次，它考虑了居住在人口少于 5 万人的聚居地的大量人口，这些人不应默认被视为农村地区。5 万人的下限在某种程度上是任意设定的，对较小聚居地的可达性在许多地区都非常重要。不同规模聚居地在资源提供、服务和机会方面的差异意味着，对于区域规划者和服务提供者而言，在考虑投资和政策对教育和健康服务、市场及就业机会获取水平影响时，对不同规模聚居地间获取情况的更细致评估将变得尤为重要。
在 2015 年，我们开发了一套更广泛的全球可访问性指标，以代表不同人口规模聚居地可获取的各种资源、服务和机会。我们使用公开可用的信息对人类聚居地进行了九级分层 11 ，从五千人或以上的小镇到五百万或以上的超大城市。对于九个级别的每一个，我们使用 gdistance R 包 13,14 来计算世界上任何地点到最近聚居地的旅行时间，具有高空间分辨率。我们使用了现有估计的穿越地球表面每 30 弧秒像素所需旅行时间的数据，使用陆地或水上最可能的交通方式 9 。输出包括九个单独的可访问性数据层，每个对应一个聚居地大小类别，表示到最近聚居地的旅行时间。
对可访问性数据层的验证是在 1,511 个 2°×2°的瓷砖上进行的，这些瓷砖代表了不同的经济资源环境，通过将我们方法论中的人类聚居点之间的旅行时间估计与 Google Maps 驾驶方向应用中的驾驶时间估计进行比较。
这些可访问性数据层提供了可持续农村发展所需的关键信息。这些层可以用于识别存在访问不平等的地区，识别可以改善访问的机会，并评估由此产生的效益。特别相关的是能够区分不同规模聚集地的旅行时间，因为越来越多的证据表明，农村人口与不同规模城市和城镇的接近度对这些农村地区的经济和社会发展的影响。
方法
旅行时间的计算基于在常规网格上移动时累积的旅行时间，网格中的每个像素都与之关联成本或旅行时间。这个网格被称为摩擦或成本表面。使用最短路径算法来找到像素与任何给定聚落大小类的附近聚落（或目标）之间的最短旅行时间，并将该最短旅行时间记录在像素中。对网格中的每个像素重复最短成本计算，以生成到最近聚落的完整旅行时间地图。为了简化处理和提高效率，地球被分割成瓷砖，最终拼接成一张地图。这个过程对每个聚落大小类重复进行，以生成九张全球地图。
输入数据集
人类人口聚居数据
我们使用了联合研究中心全球人类定居层（GHSL） 11 的 2016 年版本数据集，这些数据集以栅格格式表示低密度城市群（LDC）和高密度城市中心（HDC）（表 1）。LDC 包括城镇、郊区或小型城市区域，定义为连续单元格，人口密度至少为每平方公里 300 人 2 ，且人口至少为 5000 人。HDC 通常指的是城市或大型城市区域，定义为连续单元格，人口密度至少为每平方公里 1500 人 2 ，或者基础设施建设密度超过 50%，且人口至少为 5 万人。GHSL 数据以 GeoTIFF 栅格格式提供，分辨率为 1 公里，在世界摩尔韦德投影（坐标参考系统 EPSG:54009）下。对于 2015 年人口估计 12 （表 1），在 WGS84 30 弧秒栅格（坐标参考系统 EPSG:4326）上进行了投影，并在 LDC 和 HDC 的范围内进行求和，以确定每个定居点的 2015 年人口。然后将这些定居区域转换为多边形，其中包含定居点 ID 和人口作为属性，并投影到 WGS84。
表 1 输入数据集。
有几个基于人口的结算层级系统。我们在这里没有依赖任何一种系统，而是从 LDC 和 HDC 数据集的特征以及常用城市规模类别中推导出每种结算类别的最小人口门槛（ 17 表 2）。最大的最小门槛（类别 1）为五百万，作为对大都市（通常定义为一千万或更多）的慷慨定义，以确保全球所有地区都有此类城市。一百万（类别 2）和五十万（类别 3）由联合国使用（ 17 ）。20 万（类别 4）和 10 万（类别 5）是在 50 万（类别 6）和 5 万（类别 7）之间选择的逻辑步骤，类别 6 是 HDC 数据集中的最小结算人口，也是两个全球可达性地图中使用的相同门槛（ 9,10 ）。同样，2 万（类别 7）和 1 万（类别 8）是在 5 万（类别 9）和 5 千（类别 10）之间选择的逻辑步骤，类别 9 是 LDC 数据集中的最小结算人口，也是 Weiss 等人（ 9 ）简单认为的农村地区中居住的超过 13 亿人口的门槛。还有 10 亿人生活在 5 万到 20 万人口的城市中，无论好坏，这与生活在百万甚至更多人口的城市相比，确实有很大的不同。
表 2 结算类別、人口门槛和特性。
不采用层级系统的原因是没有单一的首选层级来定义访问权限；层级系统需要根据所研究的问题进行定制。例如，关于从农村地区日常通勤的研究可能会决定依赖于到 5 万或以上人口城市的时间阈值以下的距离来作为访问依据。从另一个角度来看，例如关于紧急医疗服务访问选项的研究，优先级可能由旅行时间与城市中心规模之间的权衡驱动，作为衡量医疗质量的代理。这里提供的九个可达性数据层将允许用户根据旅行时间和城市定居点中感兴趣的特征来自行选择如何优先考虑访问。
旅行时间数据
估算地球上任意位置到最近聚居地的旅行时间需要一个成本表面，该表面估计穿过地表上每个像素所需的时间。基于栅格的模型允许在矢量基于网络模型中不可行的非网络区域内的移动。使用了 2015 年的全球成本（摩擦）表面（分辨率 30 弧秒 9 ），其分辨率为 30 秒，用于本次分析（表 1）。像素值为每米分钟，数据格式为 GeoTIFF。
摩擦表面整合了关于运输网络（公路、铁路、河流、运河和海上航线）以及它们上的旅行速度的最新信息，以及离网区域（土地覆盖、坡度和海拔）的特征和它们上的典型步行速度，以及穿越国界所需的时间。摩擦表面将所有这些信息合并到一个层中，其中结果像素值是从所有输入中获得的最快旅行速度（因此，包含道路和河流的像素将假设道路跨越河流，像素值将代表通过道路穿越该像素所需的时间）。生成摩擦表面的过程在 Weiss 等人 9 中详细描述。
土地遮罩
全球范围内，在 85°N 和 60°S 之间，对所有陆地和海平面像素进行了旅行时间的计算。最终的旅行时间表面被屏蔽，仅显示陆地和内陆水域区域。从全球行政边界数据集（GADM）v3.6 18 的 Shapefile 版本中，获得了 30 弧秒的屏蔽。
旅行时间的计算
旅行时间计算需要表示旅程起点和终点的点位置。起点位置是成本表面中每个像素的中心。终点位置是人类聚落多边形边界上的均匀分布点。因此，每个人类聚落都有多个表示其位置的点，实际上在模型中是一个多边形，这比单个点更准确地表示了城市区域。我们的计算表示从给定位置到聚落边界上最近点的旅行时间。另一种选择，即计算到每个聚落中心的旅行时间，看起来很有吸引力，但聚落中心在 GHSL 数据集中未定义，生成中心意味着对聚落中心的任意定义，而凹多边形的中心可能位于其边界之外。
处理区域
所有旅行时间处理都在 R 中使用 gdistance 包完成 13 . 由于 gdistance 使用无法在大型栅格上操作的函数，因此处理是在分区的基础上进行的。处理的最大分区大小估计为大约 40°×40°度，使用 30 弧秒栅格。这些分区需要相邻分区之间有相当大的重叠，最多 20°，以避免在计算像素与最近目标之间的最短旅行时间时出现伪像。在试错过程中，生成了 25 个重叠分区（图 1），大小不同，用于处理，并生成了 11 个分区，所有输入数据重新定位在子午线（180°E 或 W）上，以避免在计算穿越子午线的旅行时间时出现任何不连续性。
处理步骤
所有处理都在 R 19 （版本 3.3.3）中完成，使用 RStudio（版本 1.0.143）在运行 Windows Server 2012 R2 标准版的 32 核服务器上，内存为 512Gb。处理依赖于 raster 20 （版本 2.6-7）, gdistance13（版本 1.2-2）和 rgdal 21 （版本 1.2-18）包。为了生成每个九个可达性图层，所有处理都在 30 弧秒分辨率的全球范围的栅格上进行。
1.剪裁摩擦表面和人类聚居地（目标）到每个区域的空间范围。
2.对于每个区域：
a.基于摩擦表面中每个像素的八方向旅行，生成转换矩阵。转换值通过相邻像素之间的图形连接表示。转换值可以计算为相邻像素值的平均值。然而，gdistance 预期的是传导值而不是电阻值。在这种情况下，像素值（穿越一米所需的时间）是一个电阻值。传导值计算为 1.0/电阻值。可以通过在四个（正交或冯·诺伊曼邻域）、八个（四个正交和四个对角或摩尔邻域）或 16 个方向（以国际象棋中的国王和骑士移动的组合可视化）中进行旅行来生成转换矩阵。我们选择了八个，因为它在 GIS 软件中是最常见的网格连接，并且它在 2015 年的数据集 9 中被使用。生成转换矩阵是整个过程中最消耗内存和最慢的部分。为了减少处理时间，为每个区域生成的转换矩阵保存到文件中，然后可以用于所有九个可达性层。
b.修正转换矩阵以考虑地图失真，以及网格单元之间的对角连接。转换矩阵考虑了一个像素八个可能的移动方向，但是，对角邻居之间的距离比正交邻居更远，这需要进行修正。在地理坐标系统（例如 EPSG:4326）中工作时，还需要进行另一项修正，因为东西向的连接在极地附近变得更短，因为子午线彼此接近。这两种失真都可以通过将每个传导矩阵值除以像素中心的大圆距离来纠正。
c.计算从区域内的任意像素到区域内的最近目标的累计成本，并将其保存为栅格。
图 1 重叠处理区域。
地图显示了用于处理的 25 个重叠区域。另外还生成了 11 个跨越子午线的区域（未显示）。国家边界来自 GADM v3.6。
3.将累积成本栅格拼接成单个全球栅格，使用最小函数，其中包含栅格之间的重叠像素。区域之间的大量重叠确保我们正确地将像素匹配到其最近的目标。
4.剪裁全球栅格到陆地掩模，并以 30 弧秒分辨率保存输出为整数 GeoTIFF。
数据记录
九个可访问性层在 figshare 存储库中可用，具有关于格式、时间分辨率和空间范围的适当元数据 22 。每个可访问性层都是以 WGS84（坐标参考系统 EPSG:4326）投影的 30 弧秒分辨率栅格，边界框为 85°N，180°E，60°S 和 180°W。格式为单波段 GeoTIFF，16 位无符号整数，无数据值为 65,535。像素值表示从该像素到最近的聚居地的时间（表 3）。报告了所有被分类为陆地或内陆水域区域的像素的旅行时间。图 2 显示了九个层中的三个到最近人类聚居地的全球旅行时间分钟地图。
表格 3 输出数据集。
图 2 全球可达性图层。
地图展示了 2015 年到达最近人类聚居地的旅行时间，对于九个可达性图层中的三个，（a）对于人口类别 1（5,000,000 及以上至 50,000,000 以下），（b）对于人口类别 5（100,000 及以上至 200,000 以下），以及（c）对于人口类别 9（5,000 及以上至 10,000 以下）。国家边界来自 GADM v3.6。
技术验证
我们在空间上验证了摩擦表面代表旅行时间的能力。摩擦表面是所有九个可达性层中旅行时间估计的基础，我们的预期是，摩擦层中的任何偏差都会对每个可达性层产生相同的影响。验证基于全球农村-城市映射项目（GRUMP）版本 1.01 中位于运输网络上的人类聚居点之间的旅行时间。我们计算了通过摩擦表面累计旅行时间估计的聚居点之间的旅行时间，并将其与通过 Google Maps 平台距离矩阵 API（应用程序编程接口-https://developers.google.com/maps/documentation/distance-matrix）报告的相同聚居点之间的旅程时间进行了比较。
谷歌地图平台被选作验证工具，因为它是一种广泛接受的出行时间确定方法。与我们基于栅格的方法不同，谷歌地图使用基于向量的评估，通过道路网络（每条道路都有定义的属性）来评估移动。因此，我们的结果和谷歌地图的结果都不代表“测量”（即真相）数据集。还应指出，我们的方法允许非基于道路的旅行，而谷歌地图的结果则不包括。在大多数情况下，这种差异并不重要，因为我们的模型会优先选择道路路线，因为这是最快的交通方式。然而，在偏远地区以及在计算岛屿上聚居点到像素的旅行时间时，我们的方法可能会涉及在道路外的更自由的移动（例如，谷歌地图方法会包括前往港口等待渡轮，而我们的模型则不会）。我们没有考虑由于一天中的时间或一周中的日期而导致的旅程时间变化，尽管 API 允许在提供特定出发或到达时间时进行这种考虑。
我们使用 GRUMP 而非 GHSL，因为我们的验证方法估算已知位置之间（即地名数据库中表示的点）的旅行时间。GRUMP 将定居点表示为点，而 GHSL 将定居点表示为多边形。API 需要输入位于或非常接近道路的位置，而在我们的情况下，它们还需要位于定居点内。这些输入位置本可以来自 GHSL，例如通过选择 GHSL 多边形边界上的一个点或生成中心点来代表定居点，但这种方法存在缺点。边界点可能不靠近道路，这时 API 将返回无数据。中心点不一定代表定居点的中心，而沿海地区等呈凸形的定居点可能会导致中心点位于定居点多边形之外。基于现有数据库的 GRUMP 点是一个更简单的替代方案，无需进一步的数据操作，并且它们提供了一组有效的地点，用于计算定居点之间的旅行时间。尽管 GRUMP 数据集中的聚落人口和范围随着时间的推移可能会发生变化，但它们的位置保持不变，人口和范围对我们验证方法并不相关。
我们在全球范围内生成了 2° × 2°的无重叠瓷砖，覆盖了摩擦表面的范围。我们识别了 1,511 个具有足够人口定居点（n）的瓷砖，以允许计算至少 10 个定居点之间的配对旅行时间[配对 = (n × (n-1))/2]。实际上这意味着至少有五个定居点。对于这些瓷砖中的每一个，如果有 5 到 10 个定居点，则我们选择了所有定居点。如果有超过 10 个定居点，则 i 被计算为定居点总数除以 10，然后我们选择了每 i th 个定居点，直到达到 10。使用累计成本函数和 Google API（通过 mapsapi R 包版本 0.4.0 24 ）计算了每个点对之间的旅行时间。由于假设旅行时间在两个方向上相同，因此得到的旅行时间矩阵是对称的，我们仅从矩阵中提取了下三角。在有 10 个定居点的情况下，最多可以得到 45 对，但瓷砖之间的配对数量会有所不同。在 1,511 块瓷砖上，我们计算了 47,812 次旅程中我们估计的旅行时间和 Google API 提供的旅程时间之间的差异和百分比差异。这提供了摩擦表面各位置之间累积旅行时间的空间评估。
我们的估计行程时间通常比谷歌 API 提供的要短。在各个瓷砖中，我们估计的中位行程时间为 88 分钟，介于 48 到 143 分钟的四分位数范围内，而谷歌 API 估计的中位行程时间为 106 分钟，介于 61 到 167 分钟的四分位数范围内。在所有瓷砖中，差异倾向于左侧，我们的时间估计比谷歌 API 报告的短 72%的瓷砖。中位数差异为-13.7 分钟，介于-35.5 到 2.0 分钟的四分位数范围内（图 3，面板 a），而 60%的瓷砖的绝对差异为 30 分钟或更少，80%的瓷砖的绝对差异为 60 分钟或更少（图 3，面板 c）。中位数百分比差异为-16.9%，介于-30.6%到 2.7%的四分位数范围内，而绝对百分比差异在 43%的瓷砖中为 20%或更少，在 80%的瓷砖中为 40%或更少（图 3，面板 d）。
图 3 验证输出。
(a) 估计时间差（我们的估计 - Google 行程时间），(b) 估计时间百分比差（100 × [我们的估计 - Google 行程时间]/Google 行程时间），(c) 时间差的直方图，(d) 百分比差的直方图。面板 a 和 c，面板 b 和 d 具有相同的配色方案。国家边界来自 GADM v3.6。
Full size image 验证中的空间信息显示，最小的差异出现在北美、阿根廷、欧洲、西俄罗斯、西亚、南非、澳大利亚和新西兰。这些地区的道路网络密集，速度属性很可能在谷歌数据库和我们的摩擦表面中都有良好的表示。我们旅行时间估计比谷歌的短的最大负差异发生在不同的地区，可能有不同原因。在安第斯地区，这可能与摩擦表面中通过山口的旅行速度过于乐观有关。坡度被用于惩罚摩擦表面中的步行速度，但不惩罚运输网络中的旅行速度。此外，将道路矢量化为 30 弧秒分辨率的栅格减少了道路的曲率，导致了较短的旅行时间估计。
一些沿海地区和岛屿群岛在两种旅行时间估计之间存在较差的匹配（日本、印度尼西亚和菲律宾。这些是无法在谷歌 API 中计算不同陆地之间定居点之间的驾驶行程时间的区域，导致某些点对的旅行时间存在较大的负差异。）
南亚，特别是印度和孟加拉国，都显示出了负差异，山区和沿海地区差异最大。这里的负差异可能再次归因于栅格化过程，其中密集但不一定高效的交通网络会导致栅格表示中的连接性异常高，与谷歌 API 基于矢量的表示相比。一些地区，如撒哈拉以南非洲，显示出负差异和正差异的混合模式，这表明交通网络质量或交通网络信息可用性在空间上的变异性。在中国观察到了大量正差异的集中。这是由于中国在摩擦表面缺乏交通网络信息所导致的已知问题，这意味着谷歌的旅程时间比通过摩擦表面估计的要快。
我们还观察到基于每块瓷砖的两个旅行估计的线性模型之间的良好匹配，两个估计之间的中位调整 R 2 值为 0.89，介于 0.77 到 0.96 之间，中位根均方误差为 15.8 分钟，介于 7.7 分钟到 31.4 分钟之间。整体良好的匹配以及两个旅行估计之间的一般低差异表明，摩擦表面代表了映射可达性时合理的旅行时间。
验证结果与可访问性图层一起包含在 22 中。
使用说明
无障碍层可以在许多地理信息系统或遥感软件中进行可视化和分析，例如 QGIS、GRASS、ENVI、ERDAS 或 ArcMap，以及统计和建模包如 R 或 MATLAB 中。它们也可以在基于云的工具中用于地理空间分析，如 Google Earth Engine。
九层代表到达不同人口范围的人类聚居地的旅行时间。可以通过记录各层中的最小像素值，将两层或多层合并为一层。例如，生成一个到 5,000 至 50,000 人口最近聚居地的旅行时间地图，可以通过取代表人口在 5,000 至 10,000，10,000 至 20,000 和 20,000 至 50,000 人之间的三层中的最小值来实现。
无障碍层还允许用户定义超越在层之间计算最小像素值的层次结构。当所有类别的并集等于全球人口，且任何两个类别的交集为空时，可以生成用户定义的完整层次结构。除此之外，所有逻辑一致性与实际问题相关的其他内容都由用户自行决定。
可达性层是相对于给定位置到最近目标的易达度的相对度量。虽然验证表明它们确实与典型的旅程时间相匹配，但它们不能被用来代表实际的旅行时间。摩擦表面的错误将在累积成本函数中累积，很可能，距离目标更远的位置的旅行时间与合理预期的旅行时间之间的差异将大于距离目标更近的位置。在提及较大城市到目的地的旅行时间时，如果感兴趣的位置极其偏远，应小心处理，尽管它们仍然是相对可达性合理表示。此外，该模型的一个关键假设是所有旅程都将使用最快的方式，并采取最短路径。
代码可用性
生成一系列可达性图层的 R 脚本以及执行验证的脚本可在 figshare 存储库 22 中获取。
