*These are updates. Original version from Corey Sparks, Ph.D. can be found here.

This lab complements the Lab 2 exercise using QGIS.

Here, we use tidycensus to read some tract data, learn its projection information, transform it to a new coordinate system, and measure some distance between features.

Call Packages from Library

library(tidycensus)
library(sf)
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(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, union

Note: you can run quietly=T after any library command to remove the messages in the white boxes with ## in your output, e.g., library(sf, quietly=T)

Read in Bexar county tracts

(Seem familiar? We did this chunk below before in Lab 1/Homework 2.)

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
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: Julia Kay Wolf, Ph.D. \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.

Re-Project Map into “South Central Texas” Projection

(Remember 2278 from the QGIS portion of Lab 2?)

Find other coordinate references here.

new_sa<-st_transform(sa_acs2, crs = 2278)
#Extract two tracts
twtr<-new_sa%>%
  filter(GEOID %in% c(48029181820, 48029110600))
# 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 x
head(tr_co)
## Simple feature collection with 2 features and 9 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2090862 ymin: 13704280 xmax: 2125261 ymax: 13758300
## Projected CRS: NAD83 / Texas South Central (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 (2090862 13758297)  48029   8305 15.2
## 2        37.8        15.2 POINT (2125261 13704278)  48029   5293 37.8

Measure Feet Apart

st_distance(tr_co)
## Units: [US_survey_foot]
##          1        2
## 1     0.00 64041.12
## 2 64041.12     0.00
64041.12/5280 #To get feet into miles
## [1] 12.129

(Remember in QGIS we got 12.13 miles?)

Using QGIS within R

This is another way to do the above task, by running a QGIS algorithm within R using the qgisprocess package.

NOTE: The qgisprocess will not install on its own with this 4.2.2 version of R.

See the vignette here for more on what this package is and some examples. See here to explain the following chunk for installing it.

install.packages("remotes")
remotes::install_github("paleolimbot/qgisprocess") #Select 1 for ALL; next there will be a pop-up box, select Yes.
library(qgisprocess) #load the package
## Using 'qgis_process' at 'C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat'.
## QGIS version: 3.28.2-Firenze
## Configuration loaded from 'C:\Users\xee291\AppData\Local/R-qgisprocess/R-qgisprocess/Cache/cache-0.0.0.9000.rds'
## Run `qgis_configure(use_cached_data = TRUE)` to reload cache and get more details.
## >>> If you need another installed QGIS version, run `qgis_configure()`;
##     see its documentation if you need to preset the path of qgis_process.
## - Using JSON for input serialization.
## - Using JSON for output serialization.
qgis_configure() #set up QGIS - find the executable
## getOption('qgisprocess.path') was not found.
## Sys.getenv('R_QGISPROCESS_PATH') was not found.
## Trying 'qgis_process' on PATH...
## 'qgis_process' is not available on PATH.
## Found 1 QGIS installation containing 'qgis_process':
##  C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat
## Trying command 'C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat'
## Success!
## QGIS version: 3.28.2-Firenze
## Saving configuration to 'C:\Users\xee291\AppData\Local/R-qgisprocess/R-qgisprocess/Cache/cache-0.0.0.9000.rds'
## Metadata of 1169 algorithms queried and stored in cache.
## Run `qgis_algorithms()` to see them.
## - Using JSON for input serialization.
## - Using JSON for output serialization.
# qgis_algorithms() lists all the available routines in QGIS
head(qgis_algorithms())
## # A tibble: 6 × 27
##   provider provi…¹ algor…² algor…³ algor…⁴ provi…⁵ can_c…⁶ depre…⁷ group has_k…⁸
##   <chr>    <chr>   <chr>   <chr>   <chr>   <chr>   <lgl>   <lgl>   <chr> <lgl>  
## 1 3d       QGIS (… 3d:tes… tessel… Tessel… 3d      TRUE    FALSE   Vect… FALSE  
## 2 gdal     GDAL    gdal:a… aspect  Aspect  gdal    TRUE    FALSE   Rast… FALSE  
## 3 gdal     GDAL    gdal:a… assign… Assign… gdal    TRUE    FALSE   Rast… FALSE  
## 4 gdal     GDAL    gdal:b… buffer… Buffer… gdal    TRUE    FALSE   Vect… FALSE  
## 5 gdal     GDAL    gdal:b… buildv… Build … gdal    TRUE    FALSE   Rast… FALSE  
## 6 gdal     GDAL    gdal:b… buildv… Build … gdal    TRUE    FALSE   Vect… FALSE  
## # … with 17 more variables: help_url <chr>, name <chr>,
## #   requires_matching_crs <lgl>, short_description <chr>, tags <list>,
## #   provider_can_be_activated <lgl>, default_raster_file_extension <chr>,
## #   default_vector_file_extension <chr>, provider_is_active <lgl>,
## #   provider_long_name <chr>, provider_name <chr>,
## #   supported_output_raster_extensions <list>,
## #   supported_output_table_extensions <list>, …

We can use grep to search for specific terms in the algorithms.

algs<-qgis_algorithms()
algs[grep(x = algs$algorithm, "distance"),"algorithm"]
## # A tibble: 24 × 1
##    algorithm                              
##    <chr>                                  
##  1 gdal:gridinversedistance               
##  2 gdal:gridinversedistancenearestneighbor
##  3 grass7:r.distance                      
##  4 grass7:r.grow.distance                 
##  5 grass7:v.distance                      
##  6 grass7:v.net.distance                  
##  7 native:extractwithindistance           
##  8 native:segmentizebymaxdistance         
##  9 qgis:distancematrix                    
## 10 qgis:distancetonearesthublinetohub     
## # … with 14 more rows
qgis_show_help("qgis:distancematrix")
## Distance matrix (qgis:distancematrix)
## 
## ----------------
## Description
## ----------------
## This algorithm creates a table containing a distance matrix, with distances between all the points in a points layer.
## 
## 
## ----------------
## Arguments
## ----------------
## 
## INPUT: Input point layer
##  Argument type:  source
##  Acceptable values:
##      - Path to a vector layer
## INPUT_FIELD: Input unique ID field
##  Argument type:  field
##  Acceptable values:
##      - The name of an existing field
##      - ; delimited list of existing field names
## TARGET: Target point layer
##  Argument type:  source
##  Acceptable values:
##      - Path to a vector layer
## TARGET_FIELD: Target unique ID field
##  Argument type:  field
##  Acceptable values:
##      - The name of an existing field
##      - ; delimited list of existing field names
## MATRIX_TYPE: Output matrix type
##  Default value:  0
##  Argument type:  enum
##  Available values:
##      - 0: Linear (N*k x 3) distance matrix
##      - 1: Standard (N x T) distance matrix
##      - 2: Summary distance matrix (mean, std. dev., min, max)
##  Acceptable values:
##      - Number of selected option, e.g. '1'
##      - Comma separated list of options, e.g. '1,3'
## NEAREST_POINTS: Use only the nearest (k) target points
##  Default value:  0
##  Argument type:  number
##  Acceptable values:
##      - A numeric value
## OUTPUT: Distance matrix
##  Argument type:  sink
##  Acceptable values:
##      - Path for new vector layer
## 
## ----------------
## Outputs
## ----------------
## 
## OUTPUT: <outputVector>
##  Distance matrix
out = qgis_run_algorithm(alg = "qgis:distancematrix",
               INPUT = tr_co[1,],
               INPUT_FIELD = "GEOID", 
               TARGET = tr_co[2,],
               TARGET_FIELD = "GEOID",
               MATRIX_TYPE = 0, 
               NEAREST_POINTS = 1)
## Using `OUTPUT = qgis_tmp_vector()`
## JSON input ----
## {
##   "inputs": {
##     "INPUT": "C:\\Users\\xee291\\AppData\\Local\\Temp\\RtmpWqGTBy\\file30d81a8870e5\\file30d818075f3b.gpkg",
##     "INPUT_FIELD": "GEOID",
##     "TARGET": "C:\\Users\\xee291\\AppData\\Local\\Temp\\RtmpWqGTBy\\file30d81a8870e5\\file30d8618768fb.gpkg",
##     "TARGET_FIELD": "GEOID",
##     "MATRIX_TYPE": 0,
##     "NEAREST_POINTS": 1,
##     "OUTPUT": "C:\\Users\\xee291\\AppData\\Local\\Temp\\RtmpWqGTBy\\file30d81a8870e5\\file30d824b5d60.gpkg"
##   }
## }
## 
## Running cmd.exe /c call \
##   "C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat" --json run \
##   "qgis:distancematrix" -
output_sf <- sf::read_sf(qgis_output(out, "OUTPUT"))
output_sf$Distance
## [1] 64041.12
64041.12/5280 #To get feet into miles
## [1] 12.129

NOTE: tr_co[1,] means select the first row from tr_co which is one of our two tracts. tr_co[2,] means select the 2nd row, which is our other tract.

(See it’s the same 12.13 mile distance?)