Main

An axiom for the twenty-first century is that the world is becoming increasingly connected. Although this is certainly true for electronic forms of communication, physical links between locations, and thus the time it takes to move between them, remain constrained by available infrastructure as well as physical and political impediments to travel. Eliminating disparities in accessibility is central to the Sustainable Development Goals (SDGs) set out by the United Nations1, which explicitly call for improved or universal access to key services, including education programs, health services, and banking and financial institutions. Cities are the epicentres of such activities3,4,5,6, and how easily people can reach urban areas directly affects whether crucial services can be obtained.

What constitutes accessibility is widely debated and precise definitions of this metric can be arbitrary. In this study, we operationalize accessibility in terms of travel time required to reach the nearest urban centre, defined as a contiguous area with 1,500 or more inhabitants per square kilometer or a majority of built-up land cover coincident with a population centre of at least 50,000 inhabitants7. We define accessibility using travel time as it is readily interpretable, can feasibly be generated at global scales, and is known to be a predictive metric in research domains including conservation8, food security6, trade9 and population health10,11. Furthermore, travel time better captures the opportunity cost of travel than Euclidean or network distance, and ultimately reflects the information humans use to inform transport decisions. The outcome of this research is a map that provides an actionable dataset that will support many research and policy needs. To demonstrate the map’s utility for global and local decision-making, we provide exploratory analyses examining relationships between accessibility and national-level income as well as economic prosperity, educational attainment, and healthcare utilization at the level of household clusters.

Our study responds to an increased need for fine-grained quantification of accessibility worldwide. The only previous assessment of global accessibility2 was for the year 2000, and marked advances in data quality and availability have since occurred. By anchoring our global accessibility map to 2015 (that is, the year of formal SDG adoption), we provide a baseline from which to track infrastructural improvements and urban expansion throughout the duration of the SDGs because accessibility is a precondition for many development targets. Although our results are useful in a variety of contexts, their potential impact centres around a more unifying aim: catalysing action to narrow gaps in opportunity by improving accessibility for remote populations and/or reducing disparities between populations with differing degrees of connectivity to cities.

To quantify travel time required to reach the nearest city via surface transport (air travel is not considered), we applied a similar methodology to that which had been used previously to produce the foremost existing global accessibility map2 to updated and expanded input data sets for 2015. These inputs consisted of gridded surfaces that quantify the geographical positions and salient attributes of roads, railways, rivers, water bodies, land cover types, topographical conditions (slope angle and elevation), and national borders. Roads are the primary driver of accessibility globally and also represent the most substantial advance from the previous accessibility mapping effort. Our roads dataset was created by merging Open Street Map (OSM) data with a distance-to-roads product derived from the Google roads database; these datasets were extracted in November 2016 and March 2016, respectively. The resulting roads dataset represents a global-scale synthesis of these two data sources, and the unparalleled data quality of this dataset led to a 4.8-fold increase in road pixels relative to the dataset used for the previous accessibility map2. Although roads built since 2000 contributed to the increased data volume, the primary driver of this increase was the inclusion of minor roads (for example, unpaved rural roads and exurban residential streets). The OSM database also provided the necessary information for assigning country- and road-type-specific speed data to road pixels. The Google roads data provided information critical for maintaining connectivity in areas where OSM coverage was sparse and/or fragmented owing to its piecemeal data collection approach. All input datasets were combined to create a global ‘friction surface’ at a resolution of approximately 1 by 1 km at the equator (that is, 30-arcsec resolution), effectively enumerating the generalized rates at which humans can move through each pixel of the world’s surface.

The Global Human Settlement Grid of high-density land cover (GHS-HDC)7 was used to represent cities for this research. This dataset consisted of 13,840 urban areas, an increase of 62.6% from the city points dataset used for the 2000 map2. We applied an algorithm that identified, for each pixel, the path with the least cost12 through the friction surface to any pixel defined as an urban area within the GHS-HDC. This approach ultimately produced a global accessibility map enumerating travel time to the closest city for all areas between 60° south and 85° north latitude (Fig. 1). We generated the friction surface and most of the accessibility map using the Google Earth Engine platform13, and by freely distributing our data and code our design supports the construction of bespoke accessibility maps for specific policy or programmatic priorities (for example, travel time to healthcare facilities, schools, employment centres, or markets). As with any global mapping effort, however, this study is subject to limitations that we describe in the Methods.

Figure 1: Global map of travel time to cities for 2015.
figure 1

The accessibility map has a spatial resolution of approximately 1 × 1 km, spans 60° south to 85° north latitude, and enumerates travel time to the city with the shortest associated journey.

PowerPoint slide

Full size image

Our accessibility map illustrates broad patterns of accessibility globally, capturing both the asymmetric distribution of cities and vast inequalities in infrastructural development. Highly accessible areas include those with abundant transport infrastructure and/or many spatially disaggregated cities, suggesting that accessibility to cities can be increased by improvements in infrastructure as well as polycentric urban development. Further exploration of accessibility relative to gridded population datasets14,15,16,17 shows that 80.7% of people (5.88 billion individuals) reside within one hour of cities, but accessibility is not equally distributed across the development spectrum. This disparity is evident when comparing accessibility for populations subdivided by World Bank classifications for income group and geographical region18 (Fig. 2), as 90.7% of people in high-income countries (concentrated in Europe and North America) live within one hour of a city compared to 50.9% of people in low-income countries (concentrated in sub-Saharan Africa). The relationship between national wealth and accessibility is more ambiguous for upper- and lower-middle income countries owing to the high population and large number of spatially diffuse urban centres found in northern India. This finding illustrates differences in accessibility that can readily be discerned from our map because it was produced using a globally consistent (and thus comparable) methodology. Although global summaries are informative, the accessibility map also supports fine-grained summarization and analysis (Extended Data Figs 1, 2, 3, Supplementary Information). Our map also provides a means for more nuanced characterizations of accessibility within rural populations than those that have been based upon commonly-used datasets such as binary classifications of urban versus rural land cover19 or national estimates of urban population percentages20 (Extended Data Fig. 4). As such, our map can be used to highlight major development gaps between predominantly urban and rural populations and provide a means of enumerating accessibility within rural populations along a continuum.

Figure 2: Global disparities in accessibility.
figure 2

a, b, Travel time of global populations grouped by World Bank income categories (a) and regions (b). Lines show populations aggregated in 10-min increments and then divided by the total population of each group. The inset charts show the percentage of the global population within each group.

PowerPoint slide

Source data

Full size image

To analyse subnational relationships between accessibility and socioeconomic, educational, and health measures, we used data collected by the Demographic and Health Surveys (DHS) program between 2000 and 2015. Our DHS database consisted of 66,768 household clusters, from 122 unique surveys spanning 52 countries, which were aggregated from nearly 1.77 million surveyed households. Although DHS surveys primarily cover low-to-middle-income countries, and thus do not fully represent all socioeconomic contexts, our results illustrate reciprocal relationships between accessibility and key metrics of human wellbeing in many geographical, political, and economic settings. We found a clear association between higher household wealth and greater accessibility to population centres (Fig. 3a). Similar patterns emerged for measures of educational attainment (Fig. 3b) and treatment of fever among children under five (Fig. 3c). Although exceptions occurred (for example, there were wealthy household clusters far from cities and vice versa), and the selected metrics were strongly collinear, the association between accessibility to cities and indicators of human wellbeing in low-to-middle-income countries was unequivocal.

Figure 3: Relating accessibility to human wellbeing.
figure 3

log-scale violin plots showing travel times (plus one minute) for household clusters in relation to metrics of wealth (a, n = 47,761), educational attainment (b, n = 59,686), and healthcare utilization (c, n = 39,014). The overlaid box plot hinges and colour-coding indicate the data quartiles, whiskers extend to the range of the data, violin shapes depict the data distributions, and medians and confidence intervals (CI = median ± 1.58(interquartile range/n0.5)) are displayed at the box plot centres.

PowerPoint slide

Full size image

The 2015 global accessibility map provides a high level of detail while also characterizing spatial heterogeneity in accessibility at a range of spatial scales. Our map is likely to serve as a critical input for future geospatial modelling endeavours, including those that highlight positive aspects of low accessibility, such as the protective effect that remoteness provides to wilderness areas21,22, or reinforce the need for strategic road building that avoids unnecessary environmental damage23,24. We illustrate such use through a cursory case study relating accessibility with forest loss between 2000 and 2015 (Extended Data Figs 1, 2, 3) based upon the Global Forest Change dataset25. This study shows the potential of our map for contributing to natural science research, conservation efforts, and formulation of environmental policy.

While results from our exploratory analyses do not causally link accessibility to metrics of development (for example, they cannot be used to determine whether places are more affluent because of greater accessibility or vice versa), they nevertheless illustrate the relationship between travel time and socioeconomic outcomes encompassed within the Sustainable Development Goals. Many view reducing inequalities in the accessibility of the services, institutions, and economic opportunities offered by cities as a vital pathway to sustainable development and improved livelihoods for all populations. Our analysis supports this perspective and future studies should track and evaluate the multifaceted effects that result from improved accessibility.

Methods

To model the time required for individuals to reach their most accessible city, we first quantified the speed at which humans move through the landscape. For this, we built on previous work that had integrated a number of infrastructural, political, and environmental datasets within a geographic information system (GIS)-based model2. The principle underlying this work was that all areas on Earth, represented as pixels within a 2D grid, had a cost (that is, time) associated with moving through them that we quantified as a movement speed within a cost or ‘friction’ surface. We then applied a least-cost-path algorithm12 to the friction surface in relation to a set of high-density urban points. The algorithm calculated pixel-level travel times for the optimal path between each pixel and its nearest city (that is, with the shortest journey time). From this work we ultimately produced two products: (a) an accessibility map showing travel time to urban centres, as cities are proxies for access to many goods and services that affect human wellbeing; and (b) a friction surface that underpins the accessibility map and enables the creation of custom accessibility maps from other point datasets of interest.

Accessibility mapping methodology

The datasets that we used to construct the friction surface characterize the spatial locations and properties of roads, railroads, rivers, bodies of water, topographical conditions (elevation and slope angle), land cover, and national borders. The datasets were converted into aligning grids with a 30 arcsec resolution, with the pixel values representing speeds of movement. The layers were then combined following the approach defined within an earlier accessibility mapping project2 such that the fastest mode of transport took precedence. The only exception to this logic was national borders, for which a crossing-time penalty was superimposed with priority over all other layers. The borders dataset was created from a UN global administrative units layer (GAUL) such that each border segment had a unique numerical identifier. This approach supports setting border-specific crossing times via a lookup table, however usable data do not presently exist for universally defining this parameter. As such, we used a static, one hour crossing-time penalty for all borders other than those within the Schengen and the UK–Ireland common security zones. Note that for readability the travel speeds for other input layers are provided in km h-1, but the actual units within the friction surface raster are minutes required to travel one metre.

With the exception of the rivers input, each of the datasets we used in this project and the methods we used to pre-process the data have improved considerably relative to those of the 2000 accessibility map that was produced in 20082. Two road datasets were combined for this research. The first road input layer consisted of vector data extracted from the OSM database, which was created by a user community dedicated to producing open-source, geocoded datasets of infrastructural resources. The OSM dataset was converted into a grid matching the geographic resolution and extent of the eventual friction surface. In cases for which vector features of more than one road type were present within a single pixel, the road type with the highest associated travel speed took precedence. This rasterization procedure resulted in an integer grid in which pixel values corresponded to a single road type that was subsequently linked to a speed via a lookup table, which we also derived from the OSM database. The lookup table contains the country-specific mean travel speeds associated with each available road type, as derived from attributes linked to individual roads by the OSM user community. We used this lookup table approach rather than direct assignment of road speeds because such speed-of-travel information was infrequently assigned to road vectors within the OSM database. This limitation also necessitated the creation of a global default lookup table, which we created using mean values for each road type from all countries. We applied values from the default table in cases where a country had no speed limit records for one or more road types found within it.

The second, equally important source of road data was the Google distance to roads surface. This Google dataset was also global in extent, although China and the Korean Peninsula were omitted owing to data distribution limitations. To combine the two road datasets the Google distance to roads raster was first restricted to include only pixels with values of 500 m or less, thereby approximating the 1 × 1 km rasterization of the OSM road vectors. Unlike the OSM data the resulting Google roads raster lacked road-type information. As such the OSM road-type designation took precedence if both layers contained road information for a single pixel. Where only Google road data were available, the pixels were given the default integer value corresponding to the generic ‘road’ class from OSM. When creating the friction surface, all pixels from the combined roads raster were assigned the road travel speeds from the OSM-based lookup tables. For the lookup procedure, we also used a grid of administrative units to determine each pixel’s country association.

The railroad input layer was also created from the rasterized OSM surface. Unlike the OSM roads data, however, the railroads were not differentiated by type within OSM and thus consisted of a single class with a uniform movement speed. The railroad speed used in this project was 24.3 km h−1, which was the mean value assigned to railroad vectors extracted from the OSM database.

Three datasets were used to account for travel time by water within the friction surface. River travel time was added via a global set of navigable rivers rasterized from the CIA World Data Bank II vector rivers dataset26, which was the only holdover input variable from the circa 2000 accessibility mapping endeavour2. Other options available for characterizing major rivers were explored, most notably the HydroSHEDS27 river network and Vector Map Level 0 (VMAP0)28 datasets, but we ultimately concluded that reusing the original data was warranted as neither of the alternatives was discernibly better given their associated limitations and lack of detail about which rivers were navigable. For inland water bodies, we used a newly created global surface-water occurrence dataset29, which we first aggregated from its native 30-m resolution to create a layer that enumerated the fraction of each pixel’s area that was covered by water at the resolution of the friction surface. In this procedure, all 30-m pixels within the resulting fractional surface-water dataset that were classified as water at least 80% of the time were considered permanent water, as 80% was the lowest occurrence value that we observed within ocean pixels when screening the data. The resulting fractional surface-water layer was then converted into a binary surface in which only pixels that were completely covered by permanent water were coded as a body of water amenable to be crossed by boat. The final dataset relating to water was a land–sea mask, which was used to identify ocean pixels. The movement speeds assigned to the water types within the friction surface were 10 km h−1 for rivers and lakes and 19 km h−1 for oceans. The value for rivers was based on inland travel speeds reported in the UK, Ireland, and Australia30. The ocean value was the average speed obtained from over 142 million observations of ocean-going passenger ships collected from the Automatic Identification System (AIS) and the Voluntary Observing Ship (VOS) program30.

For all pixels not covered by any of the water, road or railroad datasets, we derived a baseline speed of movement overland (that is, on foot) using the MODIS MCD12Q1 land cover product31 in which we assigned each land cover type a travel speed from a lookup table. The lookup table was created by summarizing results from an online survey designed to crowd-source estimates of how long it takes individuals to traverse each land cover type. The survey consisted of representative photos and global maps of each land cover type. Respondents were asked to estimate the amount of time it would take them to travel one kilometre (or one mile) on foot through each land cover type. The survey received 407 complete responses and, after standardizing the distance units, the median values for the fifteen land cover classes within the survey (in units of km h−1) were as follows: evergreen needleleaf forest = 3.24, evergreen broadleaf forest = 1.62, deciduous needleleaf forest = 3.24, deciduous broadleaf forest = 4.00, mixed forest = 3.24, closed shrublands = 3.00, open shrublands = 4.20, woody savannas = 4.86, savannas = 4.86, grasslands = 4.86, permanent wetlands = 2.00, croplands = 2.50, cropland/natural vegetation = 3.24, snow and ice = 1.62, and barren or sparsely vegetated = 3.00. The two land cover classes we excluded from the survey were (a) urban and built-up, which was given a speed of 5 km h−1, but this value is almost never needed owing to the higher speed (and thus precedence) of roads that dominate urban landscapes at 1 × 1 km resolution, and (b) open water, which was given a speed of 1 km h−1. The speed for the open water pixels was assigned using the rationale that if these pixels were not considered inland water within the water bodies layer (and would thus have received the inland water speed associated with boat travel) they were probably more akin to permanent wetland pixels that had a very high subpixel fraction of water to circumnavigate on foot. As such, all pixels that were classified as open water in the land cover layer but not as permanent water in the water bodies layer were given a speed half as fast as the crowd-sourced median speed for permanent wetlands of 2 km h−1.

The land-cover-dependent travel speeds were then adjusted to take into account the effect of topographical properties. Topographical data-sets used in this analysis were produced from the Global Multi-resolution Terrain Elevation Dataset 2010 (GMTED2010), a derivative of the Shuttle Radar Topography Mission data and produced by USGS32. The adjustment that we applied to elevation accounts for decreasing atmospheric density (and thus available oxygen) with altitude, which closely parallels the drop in maximal oxygen consumption (that is, VO2 max, a measure of optimal heart and lung function) and thus decreased the predicted travel speed as a function of altitude33. On the basis of the standard atmosphere calculation, equation (1) shows the adjustment factor that we associated with elevation (in metres). We treated slope angle (in degrees) similarly, as steep terrain slows humans’ ability to traverse it on foot. For the slope adjustment, we used Tobler’s Hiking Function34 as shown in equations (2) and (3), with Tobler’s walking speed capped to a maximum of 5 km h−1 and then divided by five to convert it into a fraction of maximum travel speed. The elevation and slope adjustment factors were subsequently multiplied by the land-cover-dependent travel speeds, thus lowering the speed of travel on foot and increasing the time required to traverse each associated pixel within the friction surface.

The final input for the accessibility map was the dataset of urban land cover, which was created using a layer from the Global Human Settlement (GHS) project7. This dataset was produced using a combination of satellite imagery and census data to map the spatial distribution of urban areas across the globe. To be consistent with data used in previous accessibility mapping research2, we selected the ‘high-density centres’ variant of the GHS dataset, which is defined as “contiguous cells with a density of at least 1,500 inhabitants per km2 or a density of built-up greater than 50% and a minimum of 50,000 inhabitants”. The dataset contained a total of 13,840 unique urban areas and, unlike the circa 2000 accessibility map in which cities were represented as single geographical points, cities extracted from the GHS consisted of a cluster of pixels, thus effectively representing urban areas as polygons. The switch from single-point to multiple-pixel representations of cities was operationalized by extracting each urban pixel’s centre coordinates and then applying the least-cost-path algorithm to only points on edges of urban areas. Over 400,000 high-density urban points were processed in this manner, not including points from the urban interiors, which were ignored to reduce redundant processing and later masked to have travel times of zero in the resulting accessibility map.

The friction surface was created entirely within Google Earth Engine, which was also used to create the majority of the accessibility surface. In contrast to the process used to create the friction surface, deriving the accessibility map was very computationally intensive and required a more complex processing chain. Within Earth Engine, accessibility surfaces were generated using the cumulativeCost function35, a least-cost-path function that was an experimental tool implemented specifically for this project but is now freely available within Earth Engine. By harnessing the computational power of the Google cloud-computing system the cumulativeCost function shortened the production time of the global accessibility surface from several months (when relying on local computing resources alone) to approximately two weeks. Despite reducing the production time substantially, the cumulativeCost function was still an evolving tool that was not yet capable of producing the global accessibility map in a single run or reliably producing output for latitudes above 60° if the friction surface was in geographical coordinates (that is, units of degrees latitude and longitude). As such, we created the global accessibility map by mosaicking a set of 31 tiles, 24 of which encapsulated the most computationally demanding areas and were generated within Earth Engine, and seven of which we created outside Earth Engine. The limitations of the least-cost-path function within Earth Engine at high latitudes were due to the nature of processing raster data stored in geographical coordinates because distances at high latitudes span far more degrees of longitude (and thus more pixels) than comparable distances at low latitudes. In order to parallelize computations efficiently, the Earth Engine cumulativeCost function required specification of a maximum search distance from the source points (that is, high-density urban land cover pixel centres), which we set to 1,500 km for most of the globe but reduced to 1,000 km in areas from 50° to 55° latitude owing to the afore-mentioned processing limitations at high latitudes. For latitudes above 50° we calculated accessibility tiles using the gdistance package in R36, thus ensuring an overlapping area of five-degrees latitude and providing data with which to compare the output maps from the differing sources (pixel values in these areas proved to be almost identical). We also calculated accessibility times locally for very remote islands at lower latitudes that were beyond the 1,500 km search distance threshold from their closest cities. Lastly, the cumulativeCost function in Earth Engine could not account for wrapping at ±180° longitude, so we created an alternative version of the friction surface centred at this longitude and reprocessed approximately one-fifth of the globe outside of Earth Engine to ensure that any pixels that had their closest cities on the opposite side of this ‘edge’ were ascribed accurate travel times. We then mosaicked all of the tiles together by selecting the minimum travel times for all pixels that fell within overlapping portions of multiple tiles. The result of this mosaicking operation is the global accessibility map shown in Fig. 1.

Model validation

We validated the accessibility map by comparing the travel times derived from least-cost-path calculations based on the friction surface with corresponding estimates derived from driving directions application within Google Maps (that is, comparison to travel time estimates derived using a network distance algorithm). The data source we used for validation consisted of settlement points from the Global Rural–Urban Mapping Project (GRUMP)19. Point pairs linking small settlements (that is, those with populations under 50,000 inhabitants) with their nearest city (that is, settlements with populations over 50,000 inhabitants) were processed using the friction surface approach to produce travel times akin to those within the accessibility map. After receiving special permission to automate the process of querying the Google driving directions application programming interface (API), we acquired validation travel times for each of the point pairs. A total of 53,091 validation point pairs were available after removing all coordinate pairs the API could not match. This approach limited the validation to places that fell along road networks, which precluded an assessment of the map’s accuracy in the most remote areas on Earth. However, the applied validation approach does thoroughly validate the map with respect to human populations as (a) most of the world’s people live in close proximity to a road of some variety, and (b) named points along road networks that we tested were indicative, from an accessibility perspective, of other points near roads at unnamed locations.

The validation results were encouraging, with an R2 of 0.66 and a mean absolute error (that is, the average difference between the travel times regardless of sign) of 20.7 min. The distribution of the differences between the travel times also matched our expectations as 86.5% of the point pairs had lower travel time estimates from the friction surface approach compared to the values derived from the Google API. We attributed this unequal distribution to the presence of roads within the OSM dataset that were not present within the Google data and were thus unknown to the Google API when it calculated travel time via the Google road network. Another factor that helped to explain the preponderance of lower travel times to cities derived using the friction surface was that it incorporated other forms of travel (for example, by water). This factor was particularly important for explaining point pairs with very large travel time disparities. Additional reasons why the friction surface approach tended to produce lower travel times relative to the Google API include (a) the abstraction of vector roads within raster space, which effectively shortened some roads by reducing their sinuosity; (b) the speed limit look-up table which assigned speeds to roads that may be unrealistically high (for example, if road conditions are poor); and (c) the friction surface approach assumed constant and optimal travel speeds, unlike the Google API that incorporated temporally varying delays related to traffic density (for example, rush-hour delays).

The geographical distribution of the 13.5% of point pairs for which the travel times estimated by the Google API were shorter than the corresponding times from the friction surface approach was heavily concentrated in China. This is noteworthy, because the Google roads dataset that we used to create the friction surface lacked road data for China. As such, this finding suggests that there were roads in China that the Google API used to estimate travel times that were unaccounted for within the friction surface (that is, not present within OSM). Because both the under- and overestimates were partially attributable to incomplete road network data from either OSM or Google, using a combination of these road data sources to produce the accessibility map represents a major strength of this research.

Description of map limitations

Almost all research projects that generate modelled data at global scales rely on assumptions, generalizations, and the use of best-available (even if suboptimal) datasets. An important example of this for our work is that the time it takes an individual to move through the landscape is mediated by far more factors than just infrastructure or landscape properties. Wealth, in particular, is a likely determinant of whether someone travels on foot rather than taking a vehicle and thus substantially affects accessibility on the level of the individual. As such, users are cautioned from assuming our travel time estimates are universally applicable. It should also be noted, however, that because the accessibility mapping system was developed within Earth Engine, alternative variants of the accessibility map (for example, a walking only travel time map) can readily be created.

Another caveat relates to transport by rail and water, and specifically how the least-cost-path algorithm is able to freely transition from these networks to roads or vice versa when, in reality, switching modes of transport typically requires a station or port. In our friction surface this reality is not reflected and thus the least-cost-path algorithm will occasionally utilize water and rail pixels unrealistically. Railroads are also problematic, because there is insufficient data within OSM to differentiate railroads by type and thus all railroads are assigned a relatively slow speed. As such, high-speed train travel is effectively absent from our map, although that point is largely moot when mapping accessibility to the nearest city as high-speed trains typically link large cities together (that is, to utilize such a network an individual is likely already within a city of 50,000 or more people) and are therefore similar to air travel within this context.

Including slope angle as an input layer also presents challenges because the level of detail inherent to topographical datasets depends on the spatial resolution of the elevation data used to generate such metrics. For example, data at a 1 × 1 km resolution can only reflect the slope angle at that resolution, and are likely to miss large changes in topographical relief (and thus slope angle) at finer resolutions. A related caveat is the isotropic handling of slope angle such that it will always slow down movement regardless of whether the least-cost-path is oriented uphill or downhill. The net result of these caveats is that the friction and accessibility surfaces are less reliable for off-road areas, and particularly in mountainous regions. It should also be noted that erroneous data within the global topographical dataset resulted in unrealistically high travel time estimates for a small cluster of pixels (that is, less than 50) in western Colombia.

Another known limitation of the accessibility map is that it ignores geopolitical conflicts, such as those currently occurring in Syria, where degraded infrastructure and other impediments to movement will greatly affect travel times. The relative ease with which a new friction surface can be generated using our methodology, however, would allow us to create a new friction surface and accessibility map that takes degraded infrastructure into account and thus identify areas affected by the reduced access to resources. Likewise, national borders are particularly challenging to incorporate into the accessibility map, because many borders are contested and/or unrecognized by the UN (for example, the border between Northern Cyprus and Cyprus) and thus not accounted for within the friction surface. A related challenge is borders that are effectively impermeable barriers to travel for most people (for example, the border between North and South Korea). As previously stated, there are simply no reliable data that quantify how long it takes to cross most land borders, much less the contested ones, and thus we applied a universal value that reflects the fact that most borders slow movement, particularly at road crossings, while avoiding any baseless assumptions. These factors highlight the need for better global data on border permeability and crossing times, particularly in light of ongoing policy changes related to transnational migrant flows.

Seasonal changes also present a major challenge when characterizing accessibility, particularly when they pertain to areas periodically inundated by water or covered by deep snow in which movement may be precluded for parts of the year and/or people may change their mode of transport (and thus their movement speed). Likewise, rare events such as floods and earthquakes can sever transportation links such as roads and bridges, thus markedly changing spatial accessibility patterns. Because the accessibility map was produced largely in Earth Engine, such modifications to transportation networks can be addressed by rapidly remaking the friction surface to reflect the changed reality on the ground. There are several crowd-sourced examples demonstrating how quickly such information can be collected and made available for analysis (for example, https://www.hotosm.org). A more common issue of temporal variability in accessibility pertains to public forms of transportation, which typically operate on schedules that produce delays in travel time as individuals wait for buses, trains, or ferries. Similarly, traffic congestion will slow travel times both predictably (for example, at rush hour or owing to construction) and unpredictably (for example, because of traffic accidents). As such, our accessibility should not be viewed as applicable at every moment, but rather a general estimate of accessibility during nominal travelling hours and in the absence of adverse conditions.

Exploratory analysis methodology: wealth, education, and healthcare utilization

The variables selected from the Demographic and Health Survey (DHS) Program database for exploratory analysis were household cluster measurements of the mode wealth index for heads of household, the mode educational attainment for heads of household, and the percentage of children receiving treatment for a fever (that is, healthcare utilization). The wealth and education variables were aggregated directly from questions asked within the surveys and owing to the categorical nature of these metrics, we selected the mode head of household values for analysis. The healthcare utilization metric, by contrast, was aggregated from individual-level data to provide cluster-level counts for both the numerator and denominator of the fever-treated fractions. For fever treatment this constitutes, respectively, the number of children (under five years of age) in each household cluster who received treatment for their fever divided by the total number of children within that household cluster who had a fever in the past two weeks. No statistical methods were used to predetermine sample size. Summaries of the DHS metrics relative to accessibility were depicted using violin plots (Fig. 3), which show the distribution and number of household clusters as the violin shape and area, respectively. To show the full data range these metrics were plotted using a logarithmic scale, which necessitated adding one minute to each survey cluster to plot those with travel times of zero. The added minute is reflected in Fig. 3, including the reported median and confidence interval values, the latter of which (derived as confidence intervals = median ± 1.58(interquartile range/n0.5) (ref. 37)) are quite narrow as a result of the large sample sizes.

Exploratory analysis methodology: forest loss

Despite potential beneficial aspects of short travel times to cities for humans, higher accessibility has an associated environmental cost owing to the relative ease with which humans can extract natural resources in places closer to population centres. This relationship is observable when comparing a global dataset of changes in forest density from 2000 to 201525 with travel time to cities. The initial step in this analysis was to aggregate the data for forest coverage in 2000 and forest loss from 2000 to 2015 from their native 30-m resolutions to match the 30-arcsec resolution of the accessibility surface. This process resulted in two grids quantifying the fraction of each pixel that contained forest in 2000 and experienced forest loss by 2015. By multiplying these layers by a grid of area per pixel (to convert the results to km2), binning (in 10-min increments) the pixels in both resulting layers according to their intersection with the accessibility map, and then dividing the binned totals for area of forest that experienced loss by the binned totals for area of forest in 2000, we obtained the summarized proportion of the original forest that experienced any level of loss for each travel time interval. The summarized results were then subdivided by country. Brazil and Indonesia were selected to illustrate the relationship between accessibility and forest loss using a combination of visual juxtaposition (Extended Data Figs 1, 2) and a comparison of summarizations of the national population, land area, and forest loss relative to accessibility (Extended Data Fig. 3). Our results indicate that forest loss tends to peak between 1 and 5 h from cities, or just beyond the relatively stable forest matrix associated with urban and suburban landscapes. This suggests that very close proximity to urban areas has a protective effect for forests, but a more nuanced assessment is that such areas were likely to have been heavily exploited before the year 2000 (that is, they had little readily harvestable forest remaining in 2000) and thus do not show major forest losses since the turn of the century. Geographical differences in both the magnitude and shape of the curves, however, reflect the importance of local context when interpreting results. Subnational patterns of forest loss also follow predictable patterns in both Brazil and Indonesia, with the least accessible forests showing the least amount of density loss.

Code availability

The core code used to create the accessibility map is available for download from the Malaria Atlas Project website (https://www.map.ox.ac.uk/accessibility_to_cities/).

Data availability

The accessibility map is available for visualization and/or download at http://roadlessforest.eu/map.html and https://www.map.ox.ac.uk/accessibility_to_cities/.