This notebook explores geospatial analysis with vector data in R in more detail. We go over file types, choosing a CRS, as well as an application with real world data.
Published
24 June 2024
Outline
Prerequisites
Geospatial Analysis 1
Intermediate R skills
Theoretical understanding of multiple regression
Basic geometry
Outcomes
After completing this notebook, you will be able to:
Manipulate geospatial objects in R using the sf package.
Perform geospatial operations on real world data.
Use the sf package to read in shapefile data from online sources.
Calculate distances using geospatial operations.
Understand applications of geospatial analysis in the context of economic research.
Disdier, Anne-Célia, and Keith Head. “The Puzzling Persistence of the Distance Effect on Bilateral Trade.” The Review of Economics and Statistics, 2008.
Tinbergen, J. (1962) Shaping the World Economy: Suggestions for an International Economic Policy. The Twentieth Century Fund, New York.
Introduction
In this notebook, we’ll continue our discussion of geospatial analysis with vector data. We’ll use trade data and shapefiles from Statistics Canada for our examples.
# load packageslibrary(tidyverse)
Warning: package 'tidyverse' was built under R version 4.4.3
Warning: package 'tibble' was built under R version 4.4.2
Warning: package 'tidyr' was built under R version 4.4.2
Warning: package 'readr' was built under R version 4.4.2
Warning: package 'purrr' was built under R version 4.4.3
Warning: package 'dplyr' was built under R version 4.4.3
Warning: package 'stringr' was built under R version 4.4.2
Warning: package 'forcats' was built under R version 4.4.2
Warning: package 'lubridate' was built under R version 4.4.2
Warning: package 'rmapshaper' was built under R version 4.4.3
library(crsuggest)
Warning: package 'crsuggest' was built under R version 4.4.3
Part 1: Understanding Geospatial Data
Note: geospatial files tend to be quite large and we use for loops in this notebook. It is possible that the code takes a while to run successfully.
File types and sizes
In the previous notebook, we only worked with user-generated data and with datasets embedded in the package spData. It’s likely that for your applications, you’ll have to work with external geospatial files.
The most common of geospatial file type is the ESRI Shapefile for vector data. This file has the extension .shp and typically comes in a zipped folder with at least three associated files. The sf function to read shapefiles is read_sf("filename.shp"). There are many places to find shapefiles online. Some good places to look include Geocommons, Natural Earth, and Living Atlas.
Other common file types for vector data include:
GeoJSON (.geojson): extension of JSON file, mostly used for storing latitude and longitude coordinates.
ESRI FileGDB (.gdb): stores objects created by ArcGIS (a popular software for geospatial analysis.
GeoPackage (.gpkg): lightweight, platform-independent, and compact format.
It’s important to keep in mind that these files come in varying forms and sizes. It’s very common for shapefiles to be extremely detailed and highly-defined, which can exponentially increase the file size and runtime of your analysis. The function ms_simplify() from the package rmapshaper simplifies shapefiles by reducing the number of vertices of the spatial objects while maintaining the integrity of the shapes.
Let’s read a shapefile of Canadian provinces and plot them on a map with the attribute province.
Zoom in on the province boundaries and look at how detailed they are. We’re likely not going to need this level of detail for our analysis, so let’s simplify the file to avoid large runtimes. We’ll do so by using the function ms_simplify() to reduce the number of vertices to 1% of our current number.
# run the code below to learn more about the function ms_simplify()# ?ms_simplify
# keep parameter indicates % of vertices to keepcanada_shape <-ms_simplify(shape, keep =0.01, keep_shapes =TRUE)
plot(canada_shape["province"])
Note how now we have simplified the boundaries of the provinces. We can use the function object.size() to compare how much memory we’re saving by using our canada_shape instead of the large shape object.
object.size(canada_shape)
29728 bytes
object.size(shape)
874360 bytes
Note: the object shape was already a simplified version of the high-definition shapefiles on the Statistics Canada website. Performing basic geospatial operations on those files could take multiple minutes!
Choosing an appropriate CRS
As explained in our previous notebook, a CRS is the grid of reference for our geographic objects. These grids can be either geographic, coded in spherical longitude and latitude coordinates, or projected, coded in planar (x,y) coordinates.
The manner in which we identify the CRS is through a parameter of the form AUTHORITY:CODE. In the previous notebook, we used the identifier EPSG:4326 to refer our geometric objects to the CRS WGS84, which is the most popular geographic CRS worldwide. When choosing a geographic CRS, you typically cannot go wrong with WGS84.
Choosing a projected CRS is often a tricky task - it depends on how you plan on using it. When projecting a geographic CRS, there will unavoidably be distortions to either the distances or the areas of the projected objects. The most common types of projections include:
Lambert azimuthal equal-area (LAEA) projections: preserve equal-area at all locations but distorts shapes beyond thousands of kilometers.
Azimuthal equidistant (AEQD) projections: maintain accurate straight-line distance between a point and the center point of the projection.
Lambert conformal conic (LCC) projections: the cone is set to keep distance and area properties reasonable, which is useful for regions covering thousands of kilometers.
Stereographic (STERE) projections: projections for polar regions, but with distortions on areas and distances thousands of kilometers from the center.
The package crsuggest has a function suggest_crs(), which takes a spatial object with a geographic CRS and returns a list of possible projected CRS’s that could be used for the given object.
We can use the function st_crs() to find out the CRS of our objects and st_transform() to change the CRS of the object.
# run the cell to get information from the CRSprint(st_crs(canada_shape))
Coordinate Reference System:
User input: NAD83
wkt:
GEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]]
st_crs(canada_shape)$IsGeographic
[1] TRUE
st_crs(canada_shape)$units_gdal
[1] "degree"
From the output above we can see that canada_shape is in a geographic CRS called NAD83. This geographic CRS is customized to North America and has an AUTHORITY:CODE of EPSG:4269.
Using S2 for operations in spherical coordinates
Now that we’ve looking in more detail at CRS’s, it’s time to introduce S2. S2 is a dependency of the sf package and is used for geometric operations on spherical coordinates (i.e., when the CRS is in longitude and latitude coordinates). S2 setups are especially important when calculating distances or setting buffers around objects.
This is better illustrated with an example. Let’s suppose we want to draw a buffer of distance 10 meters around the province of Alberta. Let’s perform this operation on our canada_shape dataset, first with S2 turned on and then with S2 turned off.
Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
dist is assumed to be in decimal degrees (arc_degrees).
Notice the difference? That happens because when S2 is turned off, the operation is performed in degrees and not in meters. Since the length of degrees change according to the longitude and latitude of the object, the buffers won’t be equidistant to the borders.
Now, let’s try to calculate the distance between Alberta and Ontario with S2 turned on and off.
Notice that when S2 is turned off, we get a very different number. The main takeaway is that we should use S2 when applying geospatial operations to objects with geographic coordinates.
The diagram below shows how R uses S2 for geometric operations.
Lovelace, Robin and Nowosad, Jakub and Muenchow, Jannes, (2019). Geocomputation with R.
Let’t turn S2 back on for the remainder of our analysis.
sf_use_s2(TRUE)
Part 2: Gravity Equation of Trade
Theory
In the 1660’s, Isaac Newton formalized the law of gravity: particles attract each other with a force that is proportional to the product of their masses and inversely proportional to the square of their distance. In 1962, Dutch economist Jan Tinberg adapted the law of gravity to international trade. Tinberg (1962) proposed that the aggregate trade flows between any two countries is “proportional to the gross national products of those countries and inversely proportional to the distance between them.” Mathematically, that means
\[
Trade = G \frac{Mass_{1}^\alpha Mass_{2}^\beta}{Distance^\theta}
\]
Since then, there have been numerous attempts to estimate the effects of both economic mass and distance on trade flows empirically. By taking logs of both sides of the equation, we can show that the gravity equation translates to
In this section, we’ll try to estimate the effect of distance on trade flows using data from Canadian provinces and US states.
Empirical Methods
Empirically, we can estimate the effect of distance on trade with an adaptation of the model above. For a Canadian province \(i\) and a US state \(j\):
\(T_{ij}\) is total merchandise volume between province \(i\) and state \(j\)
\(Y_{j}\) is the GDP of state \(j\)
\(X_{ij}\) is the distance between province \(i\) and state \(j\)
\(D_{i}\) is a dummy for province \(i\)
This specification is similar to the traditional gravity equation; by controlling for provinces with dummy variables, we implicitly control for the economic masses of those provinces.
The effect of interest is the coefficient on distance, \(\beta_{2}\), which can be interpreted as trade elasticity: the percentage change in total merchandise volume from a 1% change in distance.
Data Cleaning and Wrangling
We have three datasets for this analysis: canada_shape, us_states, and trade_data. The first 2 datasets contain the shapefiles of Canadian provinces and US states, and the last file contains trade data from Statistics Canada. Let’s load our datasets.
# A tibble: 6 × 5
...1 province state tmerch_value_usd gdp
<dbl> <chr> <chr> <dbl> <dbl>
1 1 Alberta Alabama 877384086 234526.
2 2 Alberta Alaska 345988655 54470.
3 3 Alberta Arizona 423062944 375545
4 4 Alberta Arkansas 121043779 132637.
5 5 Alberta California 4450662689 3062159.
6 6 Alberta Colorado 3916792877 397702.
# view canada_shape (we've loaded it in the previous section)head(canada_shape)
Simple feature collection with 6 features and 1 field
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: -141.0114 ymin: 44.99114 xmax: -57.10776 ymax: 83.10995
Geodetic CRS: NAD83
# A tibble: 6 × 2
province geometry
<chr> <GEOMETRY [°]>
1 British Columbia MULTIPOLYGON (((-128.0536 50.87683, -128.3035 50.60872, …
2 Quebec MULTIPOLYGON (((-66.32593 48.06076, -64.34198 48.41909, …
3 Nunavut MULTIPOLYGON (((-85.49845 65.92457, -86.39498 64.58634, …
4 Prince Edward Island POLYGON ((-63.98864 47.03119, -64.40965 46.70758, -63.58…
5 Saskatchewan POLYGON ((-102 60, -110 60, -110 56.92908, -110.0057 54.…
6 Yukon POLYGON ((-136.4686 68.86875, -141.0114 69.64746, -141.0…
# view us_states (comes with spData package)head(us_states)
Simple feature collection with 6 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -114.8136 ymin: 24.55868 xmax: -71.78699 ymax: 42.04964
Geodetic CRS: NAD83
GEOID NAME REGION AREA total_pop_10 total_pop_15
1 01 Alabama South 133709.27 [km^2] 4712651 4830620
2 04 Arizona West 295281.25 [km^2] 6246816 6641928
3 08 Colorado West 269573.06 [km^2] 4887061 5278906
4 09 Connecticut Norteast 12976.59 [km^2] 3545837 3593222
5 12 Florida South 151052.01 [km^2] 18511620 19645772
6 13 Georgia South 152725.21 [km^2] 9468815 10006693
geometry
1 MULTIPOLYGON (((-88.20006 3...
2 MULTIPOLYGON (((-114.7196 3...
3 MULTIPOLYGON (((-109.0501 4...
4 MULTIPOLYGON (((-73.48731 4...
5 MULTIPOLYGON (((-81.81169 2...
6 MULTIPOLYGON (((-85.60516 3...
By looking at the “trade_data” dataset, we see that:
We already have data on total merchandise volume (tmerch_value_usd) between provinces and states
We already have data on the GDP of US states (gdp)
We can easily create dummies for each Canadian province with the column province
Therefore, the only variable missing for our specification is the distance between each province-state pair. We can use our shapefiles to calculate the distance between provinces and states with the function st_distance from the sf package.
We want a dataset with the following structure:
Province
State
Total Merchandise Volume (USD)
State GDP (USD)
Distance (Km)
Alberta
Alabama
USD
USD
Km
Alberta
Arkansas
USD
USD
Km
Alberta
Arizona
USD
USD
Km
The easiest way to do this is by creating an sf object with the distance values, and then merging it with trade_data using the fields province and state.
Before we do that, let’s take a look at the CRS of our shapefiles.
# look at CRS detailsst_crs(canada_shape)$WktPretty #WktPretty makes it easier to see the CRS parameters
Our objects reference the same CRS - the geographic NAD83 CRS. We could just go ahead and calculate the distance between our objects (as long as our S2 is turned on!); however, for learning purposes, let’s take a look at the projected suggestions from the function suggest_crs().
# find list of suggested CRSsuggest_crs(canada_shape, type ="projected", units ="m") #note that the units must be in metres instead of km
suggest_crs(us_states, type ="projected", units ="m")
# A tibble: 10 × 6
crs_code crs_name crs_type crs_gcs crs_units crs_proj4
<chr> <chr> <chr> <dbl> <chr> <chr>
1 6350 NAD83(2011) / Conus Albers project… 6318 m +proj=ae…
2 5072 NAD83(NSRS2007) / Conus Albers project… 4759 m +proj=ae…
3 5071 NAD83(HARN) / Conus Albers project… 4152 m +proj=ae…
4 5070 NAD83 / Conus Albers project… 4269 m +proj=ae…
5 5069 NAD27 / Conus Albers project… 4267 m +proj=ae…
6 6500 NAD83(2011) / Minnesota Central project… 6318 m +proj=lc…
7 3594 NAD83(NSRS2007) / Minnesota Ce… project… 4759 m +proj=lc…
8 2811 NAD83(HARN) / Minnesota Central project… 4152 m +proj=lc…
9 26992 NAD83 / Minnesota Central project… 4269 m +proj=lc…
10 6504 NAD83(2011) / Minnesota South project… 6318 m +proj=lc…
This is quite an extensive list of candidates for our projected CRS. As you can see, the suggested CRS are very specific to the center of the geographic CRS - most of the suggested CRSs for Canada are artic/polar, while for the US they are based in Minnesota. Since our datasets are already in the same CRS, let’s stick with that.
# selecting the columns we'll need from `us_states`us_shape <- us_states %>%select(state = NAME, geometry)
Let’s plot our shape files on a map and then turn to merging our datasets.
Since our shapefiles are in two separate files, the approach we’re going to take is:
Create a vector with the labels province and state
Create a vector with the distance between each province-state pair
Merge those vectors together in a dataframe
Inner join the dataframe with the trade_data file
By approaching the data wrangling step-by-step, we can see exactly what is going on and can check our work each step of the way. This is often preferred when working with multiple datasets and merging through variables that are not factors (we have to merge datasets through province-state names!).
Let’s start with Step (1). We’ll create nested for loops, indexing provinces with i and states with j. Follow the comments of the code.
province_num <-nrow(canada_shape) # each entry of `canada_shape` is a provincestates_num <-nrow(us_shape) # each entry of `us_shape` is a stateprov_state_vec <-c() # create a placeholder vector that will contain the `province - state` pairsfor (i in1:province_num){ # start the loop with a province i, from a set of 1 to `province_num` provinces prov_i <- canada_shape$province[i] # select province name indexed by ifor (j in1:states_num){ # for the selected province i, loop through the states, indexed by j state_j <- us_shape$state[j] # select state name indexed by j prov_state_pair <-paste(prov_i, "-", state_j) # paste prov i and state j as a string separated by a dash prov_state_vec <-append(prov_state_vec, prov_state_pair) # append the `province-state` pair "i - j" to the placeholder vector } # close state loop for state j and repeat process for state j + 1 } # close province loop for province i and repeat process for province i + 1head(prov_state_vec)
[1] "British Columbia - Alabama" "British Columbia - Arizona"
[3] "British Columbia - Colorado" "British Columbia - Connecticut"
[5] "British Columbia - Florida" "British Columbia - Georgia"
Let’s check if our loop worked as intended. We have 13 provinces and 49 states, so we should have a total of \((13)(49) = 637\) province-state pairs.
Note: you can check the number of provinces and states by running nrow(canada_shape) and nrow(us_shape). We’re missing Hawaii and Alaska, and we’re counting the District of Columbia as a state.
length(prov_state_vec)
[1] 637
It looks like it worked! Now to Step 2. We repeat the same process, but instead of pasting strings, we calculate the distance between province i and state j. As noted earlier, our projected CRS is in meters, so those are the units of our distances.
Note: this operation might take a while to run.
province_num <-nrow(canada_shape)states_num <-nrow(us_shape)distance_vec <-c()for (i in1:province_num){ geo_prov_i <-st_make_valid(canada_shape$geometry[i]) # retrieve province geometry and adjust polygons with multiple verticesfor (j in1:states_num) { geo_state_j <-st_make_valid(us_shape$geometry[j]) # retrieve state geometry and adjust polygons with multiple vertices distance <-as.double(st_distance(geo_prov_i, geo_state_j)) # calculate distance in meters distance_vec <-append(distance_vec, distance) }}head(distance_vec)
Do you see the distances? For Step 3, let’s merge both vectors into a dataframe with cbind, and then separate the state and province match into 2 columns with separate_wider_delim.
# A tibble: 6 × 3
province state distance_m
<chr> <chr> <dbl>
1 British Columbia Alabama 2622120.
2 British Columbia Arizona 1334292.
3 British Columbia Colorado 972564.
4 British Columbia Connecticut 3213385.
5 British Columbia Florida 2987075.
6 British Columbia Georgia 2795722.
Note: The reason we are able to use cbind() to merge the two dataframes without a unique identifier is because they are ordered exactly the same. Our two for loops create pairs in an identical manner, so the distance pairs will ordered the same as our province-state pairs. We would not be able to do this is if they were ordered differently!
As a sanity check, let’s calculate the number of pairs per province (each province should be paired with exactly 49 states) and confirm that provinces and states that border each other have distance equal zero.
# group observations by `province` and calculate number of observations per provincenumber_pairs <- distance_data %>%group_by(province) %>%summarize(n =n())number_pairs
# A tibble: 13 × 2
province n
<chr> <int>
1 Alberta 49
2 British Columbia 49
3 Manitoba 49
4 New Brunswick 49
5 Newfoundland and Labrador 49
6 Northwest Territories 49
7 Nova Scotia 49
8 Nunavut 49
9 Ontario 49
10 Prince Edward Island 49
11 Quebec 49
12 Saskatchewan 49
13 Yukon 49
# find border countries for distance in meters equal zeroborder_pairs <- distance_data %>%filter(distance_m ==0)print(border_pairs)
# A tibble: 8 × 3
province state distance_m
<chr> <chr> <dbl>
1 British Columbia Washington 0
2 Quebec New York 0
3 Quebec Vermont 0
4 Quebec Maine 0
5 Quebec New Hampshire 0
6 Ontario Minnesota 0
7 Ontario New York 0
8 New Brunswick Maine 0
We only get 8 borders! That’s clearly not accurate.
What’s going on here is that whenever dealing with geospatial data, we’re always working with approximations of actual physical objects, so the data is bound to have minor errors - it’s very hard to code physical elements into just a few data points! Let’s add a tolerance of 100m to our definition of border.
Think Deeper: how would you choose a reasonable tolerance for distance calculations?
# find border countries for distance in meters equal zeroborder_pairs <- distance_data %>%filter(distance_m <100)print(border_pairs)
# A tibble: 14 × 3
province state distance_m
<chr> <chr> <dbl>
1 British Columbia Montana 20.6
2 British Columbia Washington 0
3 Quebec New York 0
4 Quebec Vermont 0
5 Quebec Maine 0
6 Quebec New Hampshire 0
7 Saskatchewan Montana 7.97
8 Saskatchewan North Dakota 26.8
9 Manitoba Minnesota 73.0
10 Manitoba North Dakota 26.8
11 Ontario Minnesota 0
12 Ontario New York 0
13 New Brunswick Maine 0
14 Alberta Montana 7.97
Much better, but we’re still missing 2 border pairs - Ontario and Michigan and British Columbia and Idaho. Let’s check their distances.
missing_pairs <- distance_data %>%filter(province %in%c("Ontario", "British Columbia"), state %in%c("Michigan", "Idaho"))missing_pairs
# A tibble: 4 × 3
province state distance_m
<chr> <chr> <dbl>
1 British Columbia Idaho 3223.
2 British Columbia Michigan 1780276.
3 Ontario Idaho 1330823.
4 Ontario Michigan 10720.
Our data says that those border pairs are a few kilometers apart. What happened here is that when we calculated the distance, we adjusted the borders of the multipolygon shapes to avoid duplicate edges with the function st_make_valid(). This affected the Ontario and British Columbia borders with some US states. Since this affects the distances in a similar manner, it likely doesn’t pose major issues to our analysis. Let’s proceed to Step 4.
# merge our `trade_data` with our `distance_data`geo_data <-merge(trade_data, distance_data, by =c("province", "state"))%>%# inner-join based on `province` and `state`mutate(dist_km =if_else(distance_m <1000, 1, distance_m/1000), # setting borders to 1 km distance to avoid log zerostmerch_mi = tmerch_value_usd/1000000,state_gdp_mi = gdp)%>%select(province, state, dist_km, state_gdp_mi, tmerch_mi)%>%filter(!is.na(dist_km), !is.na(tmerch_mi), !is.na(state_gdp_mi)) # remove pairs with no trade flowshead(geo_data)
province state dist_km state_gdp_mi tmerch_mi
1 Alberta Alabama 2364.9071 234526.4 877.3841
2 Alberta Arizona 1334.2918 375545.0 423.0629
3 Alberta Arkansas 1866.5300 132637.2 121.0438
4 Alberta California 905.1146 3062158.9 4450.6627
5 Alberta Colorado 892.6367 397701.6 3916.7929
6 Alberta Connecticut 2916.9309 285466.4 111.5382
Great! We finally compiled our data. Let’s look at how many pairs we have.
# group observations by `province` and calculate number of observations per provincenumber_pairs_final <- geo_data %>%group_by(province) %>%summarize(n =n())number_pairs_final
# A tibble: 13 × 2
province n
<chr> <int>
1 Alberta 49
2 British Columbia 49
3 Manitoba 49
4 New Brunswick 49
5 Newfoundland and Labrador 48
6 Northwest Territories 15
7 Nova Scotia 49
8 Nunavut 8
9 Ontario 49
10 Prince Edward Island 49
11 Quebec 49
12 Saskatchewan 49
13 Yukon 36
It looks like there are a few province-country pairs with few trade flows. Let’s filter out provinces that don’t trade with all 49 US states.
Think Deeper: how would not removing those pairs affect the results of our regression? Instead of filtering our dataset, what else could we do to solve this issue?
Exploring our Data
Now that we have our data, let’s plot a scatterplot to visualize our relationships. Let’s filter the data to only British Columbia pairs and plot the distance between the BC-state pairs against the tmerch_mi trade between BC-state pairs. Let’s also set the size of each data point to be the GDP of the US states.
# filter data to BC-state pairsBC_data <- model_data %>%filter(province =="British Columbia")# plot dist_km against log(tmerch_mi) and add a line of best fitggplot(BC_data, aes(x = dist_km, y =log(tmerch_mi), size = state_gdp_mi)) +geom_point() +geom_smooth(method = lm, se =FALSE, linewidth =1)
`geom_smooth()` using formula = 'y ~ x'
The line of best fit shows a negative relationship between distance and trade. Also note that the largest data points are above the line: all else equal, higher GDP is associated with higher trade levels.
Running our Model
Now, let’s run our model using lm().
model1 <-lm(log(tmerch_mi) ~log(dist_km) + province +log(state_gdp_mi), data = model_data)coeftest(model1, vcov = vcovHC)
Our regression indicates that a 1% increase in distance is associated with a 0.5% decrease in trade flows. Our parameter is statistically significant at the 1% level. As expected, the state GDP control is also positive and significant.
Discussion
Our geospatial analysis suggests that the relationships drawn by Tilberg hold for Canada and the US: distance is inversely proportional to trade flows. The negative relationship between distance and trade flows has been consistently found in the literature, but the magnitude of the effects vary with time and geography.
In a comprehensive study about the effects of distance on bilateral trade, Disdier and Head (2008) compiled 1467 different distance effect estimates from 103 papers. They found that the mean distance effect is 0.9 - that is, a 10% increase in distance lowers bilateral trade by about 9%. Furthermore, 90% of the estimates compiled fall within the range of 0.28 to 1.55, dependent on the time period and geography of the study. It seems that our results are in line with what the researchers found, which adds to the credibility of our estimates from a simple model.
This example concludes our modules on geospatial analysis. If you’re curious to learn more about the topics explored in these modules, refer to: