DEM 7093 Homework 3 using R - Map projections and transformations
Brian Surratt
10 February, 2023
Call Packages from Library
library(tidycensus)
library(sf)## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionRead in Bexar county tracts
sa_acs <- get_acs(geography = "tract",
state="TX",
county = c("Bexar"),
year = 2019,
variables=c( "DP05_0001E",
"DP03_0119PE"),
geometry = T,
output = "wide")## Getting data from the 2015-2019 5-year ACS## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## Using the ACS Data Profile# create a county FIPS code - 5 digit
# substr() - extract or replace substrings in a character vector
sa_acs$county<-substr(sa_acs$GEOID, 1, 5)
# rename variables and filter missing cases
sa_acs2 <- sa_acs %>%
mutate(totpop = DP05_0001E, ppov=DP03_0119PE) %>%
# st_transform(crs = 102740) %>%
na.omit()Find Coordinate System of Current Map
See Coordinate Systems, Projections, and Transformations from ArcGIS Pro/ESRI for more information about these concepts.
st_crs(sa_acs2)## 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]]We see these tracts are in a Geographic Coordinate System (GCS) called North American Datum 1983 (NAD83).
Create Basic Map
library(tmap)
library(tmaptools)
tm_shape(sa_acs2) +
tm_polygons("ppov",
title="% in Poverty",
palette="Blues",
style="quantile",
n=5 )+
tm_format("World",
main.title="San Antonio Poverty Estimates (2019) - Quintile Breaks",
main.title.position=c('center','top'),
main.title.size=1.5,
title="Author: Brian Surratt \nSource: ACS 2019",
legend.title.size=1.7,
legend.outside=T,
legend.text.size=1.2)+
tm_scale_bar(position = c("left","bottom"))+
tm_compass()
Click here
for more tmap aesthetic features.
Click here for a quick discussion on quantile vs. quintile.
First find the distance without doing a projection.
#Extract two tracts
twtr <- sa_acs2 %>%
filter(GEOID %in% c(48029181820, 48029110600)) #FIPS codes for the campuses
# get centroid coordinates for two tracts (these two tracts are where UTSA Main and Downtown Campuses are)
tr_co <- st_centroid(twtr)## Warning in st_centroid.sf(twtr): st_centroid assumes attributes are constant
## over geometries of xhead(tr_co)## Simple feature collection with 2 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -98.61502 ymin: 29.43017 xmax: -98.50751 ymax: 29.57909
## Geodetic CRS: NAD83
## GEOID NAME DP05_0001E DP05_0001M
## 1 48029181820 Census Tract 1818.20, Bexar County, Texas 8305 833
## 2 48029110600 Census Tract 1106, Bexar County, Texas 5293 423
## DP03_0119PE DP03_0119PM geometry county totpop ppov
## 1 15.2 9.1 POINT (-98.61502 29.57909) 48029 8305 15.2
## 2 37.8 15.2 POINT (-98.50751 29.43017) 48029 5293 37.8Measure Distance
st_distance(tr_co)## Units: [m]
## [,1] [,2]
## [1,] 0.00 19555.86
## [2,] 19555.86 0.0019555.86*3.28/5280 #To get meters into miles## [1] 12.14834(In QGIS we got 12.13 miles.)
Find other coordinate references here.
Re-Project Map into “EPSG: 3665: NAD83 / Texas Centric Albers Equal Area” projection
new_sa <- st_transform(sa_acs2, crs = 3665)Run map again to compare shape before and after projection.
tm_shape(new_sa) +
tm_polygons("ppov",
title="% in Poverty",
palette="Blues",
style="quantile",
n=5 )+
tm_format("World",
main.title="San Antonio Poverty Estimates (2019) - Quintile Breaks",
main.title.position=c('center','top'),
main.title.size=1.5,
title="Author: Brian Surratt \nSource: ACS 2019",
legend.title.size=1.7,
legend.outside=T,
legend.text.size=1.2)+
tm_scale_bar(position = c("left","bottom"))+
tm_compass()
#Extract two tracts
twtr <- new_sa %>%
filter(GEOID %in% c(48029181820, 48029110600)) #FIPS codes for the campuses
# get centroid coordinates for two tracts (these two tracts are where UTSA Main and Downtown Campuses are)
tr_co <- st_centroid(twtr)## Warning in st_centroid.sf(twtr): st_centroid assumes attributes are constant
## over geometries of xhead(tr_co)## Simple feature collection with 2 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1633963 ymin: 7257682 xmax: 1644584 ymax: 7274079
## Projected CRS: NAD83(NSRS2007) / Texas Centric Albers Equal Area
## GEOID NAME DP05_0001E DP05_0001M
## 1 48029181820 Census Tract 1818.20, Bexar County, Texas 8305 833
## 2 48029110600 Census Tract 1106, Bexar County, Texas 5293 423
## DP03_0119PE DP03_0119PM geometry county totpop ppov
## 1 15.2 9.1 POINT (1633963 7274079) 48029 8305 15.2
## 2 37.8 15.2 POINT (1644584 7257682) 48029 5293 37.8Measure Distance
st_distance(tr_co)## Units: [m]
## 1 2
## 1 0.00 19536.16
## 2 19536.16 0.0019536.16*3.28/5280 #To get feet into miles## [1] 12.1361(In QGIS we got 12.13 miles.)
Trying a different projection. This projection is SR-ORG: 7629, Bexar County, TX.
new_sa2 <- st_transform(sa_acs2, crs = 7629)#Extract two tracts
twtr2 <- new_sa2 %>%
filter(GEOID %in% c(48029181820, 48029110600)) #FIPS codes for the campuses
# get centroid coordinates for two tracts (these two tracts are where UTSA Main and Downtown Campuses are)
tr_co2 <- st_centroid(twtr2)## Warning in st_centroid.sf(twtr2): st_centroid assumes attributes are constant
## over geometries of xhead(tr_co2)## Simple feature collection with 2 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -2560830 ymin: -4433127 xmax: -2530876 ymax: -4375751
## Projected CRS: NAD83(2011) / WISCRS Rock (ftUS)
## GEOID NAME DP05_0001E DP05_0001M
## 1 48029181820 Census Tract 1818.20, Bexar County, Texas 8305 833
## 2 48029110600 Census Tract 1106, Bexar County, Texas 5293 423
## DP03_0119PE DP03_0119PM geometry county totpop ppov
## 1 15.2 9.1 POINT (-2560830 -4375751) 48029 8305 15.2
## 2 37.8 15.2 POINT (-2530876 -4433127) 48029 5293 37.8Measure Distance
st_distance(tr_co2)## Units: [US_survey_foot]
## 1 2
## 1 0.00 64724.87
## 2 64724.87 0.0064724.87/5280 #To get feet into miles## [1] 12.2585(In QGIS we got 12.13 miles.)
Trying a different projection. This projection is EPSG: 3572, North Pole LAEA Alaska.
new_sa3 <- st_transform(sa_acs2, crs = 3572)#Extract two tracts
twtr3 <- new_sa3 %>%
filter(GEOID %in% c(48029181820, 48029110600)) #FIPS codes for the campuses
# get centroid coordinates for two tracts (these two tracts are where UTSA Main and Downtown Campuses are)
tr_co3 <- st_centroid(twtr3)## Warning in st_centroid.sf(twtr3): st_centroid assumes attributes are constant
## over geometries of xhead(tr_co3)## Simple feature collection with 2 features and 9 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 5017928 ymin: -4007914 xmax: 5036601 ymax: -4007371
## Projected CRS: WGS 84 / North Pole LAEA Alaska
## GEOID NAME DP05_0001E DP05_0001M
## 1 48029181820 Census Tract 1818.20, Bexar County, Texas 8305 833
## 2 48029110600 Census Tract 1106, Bexar County, Texas 5293 423
## DP03_0119PE DP03_0119PM geometry county totpop ppov
## 1 15.2 9.1 POINT (5017928 -4007914) 48029 8305 15.2
## 2 37.8 15.2 POINT (5036601 -4007371) 48029 5293 37.8Measure Distance
st_distance(tr_co3)## Units: [m]
## 1 2
## 1 0.00 18680.77
## 2 18680.77 0.0018680.77*3.28/5280 #To get meters into miles## [1] 11.60472(In QGIS we got 12.13 miles.)
In general, why is it important to have an accurate system of projection?
Because the earth is roughly spherical, geographic coordinate systems must incorporate angular degrees. Methods for displaying geographic coordinate systems on flat surfaces are known as projections. Because of the differences in shape (sphere versus plane) mathematical calculations must be used in the conversion. There is no perfect way to display a curved surface on a flat surface without distortion. There are various ways to do it using different mathematical algorithms with different properties. So, a projection that emphasizes the accuracy of distances may compromise in the accuracy of shapes. Therefore, it is important to use a projection that preserves the feature that is most important to your project. It is important to have an accurate system of projection in order to obtain the most accurate measures when working with a flat representation of a map.
How could your results be sensitive to this?
Inaccurate projections create error in geographic measurements. GIS tools like QGIS or R are working with flat representations of data that originally describes a curved space. Due to the mathematical complexity of the process, there are multiple projection methods. Projections can be designed to preserve features of a coordinate system like distance or shape. Or projections can be designed for particular places. When applying a projection for particular features in a particular place (like measuring distances in Texas), error can be introduced by using a projection that preserves a different feature, such as shapes, or a projection designed for a particular place, like the north pole.