DEM 5093/7093 - MHW6 - apping NAICS codes using AtoZ database
Ricardo Henrique Lowe, Jr., M.S.
24 March, 2021
devtools::install_github## Error in get(genname, envir = envir) : object 'testthat_print' not found## function (repo, ref = "HEAD", subdir = NULL, auth_token = github_pat(quiet),
## host = "api.github.com", dependencies = NA, upgrade = c("default",
## "ask", "always", "never"), force = FALSE, quiet = FALSE,
## build = TRUE, build_opts = c("--no-resave-data", "--no-manual",
## "--no-build-vignettes"), build_manual = FALSE, build_vignettes = FALSE,
## repos = getOption("repos"), type = getOption("pkgType"),
## ...)
## pkgbuild::with_build_tools({
## ellipsis::check_dots_used(action = getOption("devtools.ellipsis_action",
## rlang::warn))
## {
## remotes <- lapply(repo, github_remote, ref = ref, subdir = subdir,
## auth_token = auth_token, host = host)
## install_remotes(remotes, auth_token = auth_token, host = host,
## dependencies = dependencies, upgrade = upgrade, force = force,
## quiet = quiet, build = build, build_opts = build_opts,
## build_manual = build_manual, build_vignettes = build_vignettes,
## repos = repos, type = type, ...)
## }
## }, required = FALSE)
## <bytecode: 0x7f9652dc34a8>
## <environment: namespace:remotes>remotes::install_github("r-spatial/mapview")## Skipping install of 'mapview' from a github remote, the SHA1 (b96de520) has not changed since last install.
## Use `force = TRUE` to force installationlibrary(mapview)
library(sf)## Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1library(censusxy)
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, unionlibrary(readxl)
satxgroc <- read_excel("/Users/rhl548/Desktop/Cosas/rmd/IPUMS International/satxgroc.xlsx")## New names:
## * Source -> Source...1
## * Source -> Source...18addr<-satxgroc[c(6, 7:9)]
names(addr)<-c("street", "city", "st", "zip")
head(addr)## # A tibble: 6 x 4
## street city st zip
## <chr> <chr> <chr> <dbl>
## 1 646 S Flores St San Antonio TX 78204
## 2 11122 Nacogdoches Rd San Antonio TX 78217
## 3 1314 Fredericksburg Rd San Antonio TX 78201
## 4 1545 S San Marcos San Antonio TX 78207
## 5 340 Enrique M Barrera Pkwy San Antonio TX 78237
## 6 1430 Austin Hwy San Antonio TX 78209results<-cxy_geocode(addr,
street = "street",
city = "city",
state ="st",
zip = "zip",
class="sf",
output = "simple")## 141 rows removed to create an sf object. These were addresses that the geocoder could not match.results.proj<-st_transform(results,
crs = 2278)mapview(results.proj)mean feature - average of coordinates
mean_feature<-apply(st_coordinates(results.proj), MARGIN = 2, FUN = mean)
mean_feature<-data.frame(place="meanfeature", x=mean_feature[1], y= mean_feature[2])
mean_feature<-st_as_sf(mean_feature, coords = c("x", "y"), crs= 2278)
mapview(mean_feature, col.regions="red")+mapview( results)Central feature - Median of coordinates
median_feature<-apply(st_coordinates(results.proj), MARGIN = 2, FUN = median)
median_feature<-data.frame(place="medianfeature", x=median_feature[1], y= median_feature[2])
median_feature<-st_as_sf(median_feature, coords = c("x", "y"), crs= 2278)
mapview(median_feature, col.regions="green")+
mapview(mean_feature, col.regions="red")+
mapview( results)Buffer points
grocbuff<- st_buffer(results.proj, dist = 2500)
mapview(grocbuff)+mapview(results.proj, col.regions="green")Convex hull plot
chull <- st_convex_hull(st_union(results))## although coordinates are longitude/latitude, st_union assumes that they are planarmapview(chull)+
mapview(results, col.regions = "green")kernel density - You need projected data for this to work right
R can do kernel density maps, but using simple features it’s kind of complicated. I will use Qgis through R instead using the qgisprocess package
library(qgisprocess)## Using 'qgis_process' at '/Applications/QGIS.app/Contents/MacOS/bin/qgis_process'.
## Run `qgis_configure()` for details.qgis_configure()## getOption('qgisprocess.path') was not found.## Sys.getenv('R_QGISPROCESS_PATH') was not found.## Trying 'qgis_process' on PATH## Error in rethrow_call(c_processx_exec, command, c(command, args), stdin, : cannot start processx process 'qgis_process' (system error 2, No such file or directory) @unix/processx.c:592 (processx_exec)## Found 1 QGIS installation containing 'qgis_process':
## /Applications/QGIS.app/Contents/MacOS/bin/qgis_process## Trying command '/Applications/QGIS.app/Contents/MacOS/bin/qgis_process'## Success!To use this, we need to find the name of the Qgis algorithm we want. qgis_algorithms() can return all available algorithms, then we can either filter it with View() or use grep to search for one.
algs<-qgis_algorithms()
algs[grepl(pattern = "density", x = algs$algorithm ),]## # A tibble: 8 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 grass7 GRASS grass7:r.li.ed… r.li.edgedensity r.li.edgedensity
## 2 grass7 GRASS grass7:r.li.ed… r.li.edgedensit… r.li.edgedensity.as…
## 3 grass7 GRASS grass7:r.li.pa… r.li.patchdensi… r.li.patchdensity
## 4 grass7 GRASS grass7:r.li.pa… r.li.patchdensi… r.li.patchdensity.a…
## 5 native QGIS (native c… native:lineden… linedensity Line density
## 6 qgis QGIS qgis:heatmapke… heatmapkernelde… Heatmap (Kernel Den…
## 7 saga SAGA saga:fragmenta… fragmentationcl… Fragmentation class…
## 8 saga SAGA saga:kernelden… kerneldensityes… Kernel density esti…qgis_show_help("qgis:heatmapkerneldensityestimation")## Heatmap (Kernel Density Estimation) (qgis:heatmapkerneldensityestimation)
##
## ----------------
## Description
## ----------------
##
## ----------------
## Arguments
## ----------------
##
## INPUT: Point layer
## Argument type: source
## Acceptable values:
## - Path to a vector layer
## RADIUS: Radius
## Argument type: distance
## Acceptable values:
## - A numeric value
## RADIUS_FIELD: Radius from field
## Argument type: field
## Acceptable values:
## - The name of an existing field
## - ; delimited list of existing field names
## PIXEL_SIZE: Output raster size
## Argument type: number
## Acceptable values:
## - A numeric value
## WEIGHT_FIELD: Weight from field
## Argument type: field
## Acceptable values:
## - The name of an existing field
## - ; delimited list of existing field names
## KERNEL: Kernel shape
## Argument type: enum
## Available values:
## - 0: Quartic
## - 1: Triangular
## - 2: Uniform
## - 3: Triweight
## - 4: Epanechnikov
## Acceptable values:
## - Number of selected option, e.g. '1'
## - Comma separated list of options, e.g. '1,3'
## DECAY: Decay ratio (Triangular kernels only)
## Argument type: number
## Acceptable values:
## - A numeric value
## OUTPUT_VALUE: Output value scaling
## Argument type: enum
## Available values:
## - 0: Raw
## - 1: Scaled
## Acceptable values:
## - Number of selected option, e.g. '1'
## - Comma separated list of options, e.g. '1,3'
## OUTPUT: Heatmap
## Argument type: rasterDestination
## Acceptable values:
## - Path for new raster layer
##
## ----------------
## Outputs
## ----------------
##
## OUTPUT: <outputRaster>
## HeatmapRun the algorithm
groc_dens<-qgis_run_algorithm(algorithm ="qgis:heatmapkerneldensityestimation",
INPUT=results.proj,
RADIUS = 5000,
PIXEL_SIZE = 100,
KERNEL = 0,
OUTPUT=file.path(getwd(), "grocdenst.TIF"),
load_output = TRUE)## Running /Applications/QGIS.app/Contents/MacOS/bin/qgis_process run \
## 'qgis:heatmapkerneldensityestimation' \
## '--INPUT=/var/folders/_1/tkj8dh2n73s3c5fp80h0td380000gn/T//RtmpYVQfSd/filec1ea341d06a1/filec1ea71865750.gpkg' \
## '--RADIUS=5000' '--PIXEL_SIZE=100' '--KERNEL=0' '--OUTPUT_VALUE=0' \
## '--OUTPUT=/Users/rhl548/Desktop/Cosas/rmd/screenshots/grocdenst.TIF'
##
## ----------------
## Inputs
## ----------------
##
## INPUT: /var/folders/_1/tkj8dh2n73s3c5fp80h0td380000gn/T//RtmpYVQfSd/filec1ea341d06a1/filec1ea71865750.gpkg
## KERNEL: 0
## OUTPUT: /Users/rhl548/Desktop/Cosas/rmd/screenshots/grocdenst.TIF
## OUTPUT_VALUE: 0
## PIXEL_SIZE: 100
## RADIUS: 5000
##
##
## 0...10...20...30...40...50...60...70...80...90...
## ----------------
## Results
## ----------------
##
## OUTPUT: /Users/rhl548/Desktop/Cosas/rmd/screenshots/grocdenst.TIFlibrary(raster)
library(RColorBrewer)
result<- qgis_as_raster(groc_dens)
projection(result)<-crs(results.proj)
mapview(result)+mapview(results.proj)Spatial join
A spatial join can combine attributes of one layer with another layer. Here I combine census variables with the groc clinic points.
library(tidycensus)
library(dplyr)
#load census tract data
sa_acs<-get_acs(geography = "tract",
state="TX",
county = "Bexar",
year = 2019,
variables=c( "DP05_0001E", "DP03_0009P", "DP03_0062E", "DP03_0119PE",
"DP05_0001E","DP02_0009PE","DP02_0008PE","DP02_0040E","DP02_0038E",
"DP02_0066PE","DP02_0067PE","DP02_0080PE","DP02_0092PE",
"DP03_0005PE","DP03_0028PE","DP03_0062E","DP03_0099PE","DP03_0101PE",
"DP03_0119PE","DP04_0046PE","DP05_0072PE","DP05_0073PE",
"DP05_0066PE", "DP05_0072PE", "DP02_0113PE") ,
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
## Using the ACS Data Profile#rename variables and filter missing cases
sa_acs2<-sa_acs%>%
mutate(totpop= DP05_0001E, pwhite=DP05_0072PE,
pblack=DP05_0073PE , phisp=DP05_0066PE,
phsormore=DP02_0066PE,punemp=DP03_0009PE, medhhinc=DP03_0062E,
ppov=DP03_0119PE)%>%
dplyr::select(GEOID, totpop, pblack, pwhite, phisp, punemp, medhhinc, ppov)
sa_acs2<-st_transform(sa_acs2, crs = 2278)
sa_trol<-st_cast(sa_acs2, "MULTILINESTRING")spjoin<-st_join(results.proj, sa_acs2)
head(spjoin)## Simple feature collection with 6 features and 12 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 2100083 ymin: 13676720 xmax: 2147585 ymax: 13769830
## Projected CRS: NAD83 / Texas South Central (ftUS)
## street city st zip GEOID totpop pblack
## 386 646 S Flores St San Antonio TX 78204 48029192100 2299 0.7
## 83 1430 Austin Hwy San Antonio TX 78209 48029120600 5850 1.4
## 97 1515 N Loop 1604 E San Antonio TX 78232 48029191817 12785 1.7
## 311 5025 Nw Loop 410 San Antonio TX 78229 48029180702 6892 0.8
## 457 8500 Jones Maltsberger Rd San Antonio TX 78216 48029120701 6416 3.1
## 36 1200 Se Military Dr San Antonio TX 78214 48029192200 2790 1.3
## pwhite phisp punemp medhhinc ppov geometry
## 386 42.4 0.0 3.8 73194 1.2 POINT (2128476 13699880)
## 83 24.5 0.7 6.2 54279 7.9 POINT (2147585 13726925)
## 97 34.2 0.3 5.6 98311 3.5 POINT (2136320 13769833)
## 311 45.6 0.1 8.0 37899 19.6 POINT (2100083 13725095)
## 457 47.2 1.6 7.0 51633 6.3 POINT (2132875 13735344)
## 36 60.4 1.3 5.8 41623 12.5 POINT (2134829 13676718)mapview(spjoin["punemp"])+mapview(sa_trol)Count points in polygons
Point in polygon operations are actually a spatial intersection (more on this next week!) where we see how many points fall within a given polygon.
sa_acs2$ngroc<- lengths(st_intersects(sa_acs2, results.proj))
mapview(sa_acs2, zcol="ngroc")+
mapview(results.proj, col.regions = "green")Thiessen/Voronoi Polygons
Thiessen or Voronoi polygons are a process where we can convert points into polygons.
algs[grepl(pattern = "voronoi", x = algs$algorithm ),]## # A tibble: 3 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 grass7 GRASS grass7:v.voronoi v.voronoi v.voronoi
## 2 grass7 GRASS grass7:v.voronoi.sk… v.voronoi.skele… v.voronoi.skele…
## 3 qgis QGIS qgis:voronoipolygons voronoipolygons Voronoi polygonsqgis_show_help("qgis:voronoipolygons")## Voronoi polygons (qgis:voronoipolygons)
##
## ----------------
## Description
## ----------------
## This algorithm takes a points layer and generates a polygon layer containing the voronoi polygons corresponding to those input points.
##
##
## ----------------
## Arguments
## ----------------
##
## INPUT: Input layer
## Argument type: source
## Acceptable values:
## - Path to a vector layer
## BUFFER: Buffer region (% of extent)
## Argument type: number
## Acceptable values:
## - A numeric value
## OUTPUT: Voronoi polygons
## Argument type: sink
## Acceptable values:
## - Path for new vector layer
##
## ----------------
## Outputs
## ----------------
##
## OUTPUT: <outputVector>
## Voronoi polygonsgroc_von<-qgis_run_algorithm(alg="qgis:voronoipolygons",
INPUT=results.proj,
OUTPUT=file.path(tempdir(), "grocvon.shp"),
load_output = TRUE)
groc_von<-sf::read_sf(qgis_output(groc_von, "OUTPUT"))
mapview(groc_von, alpha=.75)+
mapview(results.proj, col.regions="green")Nearest Neighbor analysis
library(spatstat)
groc.pp<-as.ppp(as(results.proj, "Spatial"))
plot(nearest.neighbour(groc.pp))
algs[grepl(pattern = "nearest", x = algs$algorithm ),]## # A tibble: 10 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 gdal GDAL gdal:gridinver… gridinversedist… Grid (IDW with nea…
## 2 gdal GDAL gdal:gridneare… gridnearestneig… Grid (Nearest neig…
## 3 native QGIS (native c… native:angleto… angletonearest Align points to fe…
## 4 native QGIS (native c… native:joinbyn… joinbynearest Join attributes by…
## 5 native QGIS (native c… native:nearest… nearestneighbou… Nearest neighbour …
## 6 qgis QGIS qgis:distancet… distancetoneare… Distance to neares…
## 7 qgis QGIS qgis:distancet… distancetoneare… Distance to neares…
## 8 qgis QGIS qgis:knearestc… knearestconcave… Concave hull (k-ne…
## 9 saga SAGA saga:knearestn… knearestneighbo… K-nearest neighbou…
## 10 saga SAGA saga:nearestne… nearestneighbour Nearest neighbourqgis_show_help("native:nearestneighbouranalysis")## Nearest neighbour analysis (native:nearestneighbouranalysis)
##
## ----------------
## Description
## ----------------
## This algorithm performs nearest neighbor analysis for a point layer.
##
## The output describes how the data are distributed (clustered, randomly or distributed).
##
## Output is generated as an HTML file with the computed statistical values.
##
## ----------------
## Arguments
## ----------------
##
## INPUT: Input layer
## Argument type: source
## Acceptable values:
## - Path to a vector layer
## OUTPUT_HTML_FILE: Nearest neighbour
## Argument type: fileDestination
## Acceptable values:
## - Path for new file
##
## ----------------
## Outputs
## ----------------
##
## OUTPUT_HTML_FILE: <outputHtml>
## Nearest neighbour
## OBSERVED_MD: <outputNumber>
## Observed mean distance
## EXPECTED_MD: <outputNumber>
## Expected mean distance
## NN_INDEX: <outputNumber>
## Nearest neighbour index
## POINT_COUNT: <outputNumber>
## Number of points
## Z_SCORE: <outputNumber>
## Z-scoregroc_nn<-qgis_run_algorithm(alg="native:nearestneighbouranalysis",
INPUT=results.proj,
OUTPUT_HTML_FILE=file.path(tempdir(), "grocnn.html"),
load_output = TRUE)## Ignoring unknown input 'load_output'## Running /Applications/QGIS.app/Contents/MacOS/bin/qgis_process run \
## 'native:nearestneighbouranalysis' \
## '--INPUT=/var/folders/_1/tkj8dh2n73s3c5fp80h0td380000gn/T//RtmpYVQfSd/filec1ea341d06a1/filec1ea2804cf5a.gpkg' \
## '--OUTPUT_HTML_FILE=/var/folders/_1/tkj8dh2n73s3c5fp80h0td380000gn/T//RtmpYVQfSd/grocnn.html'
##
## ----------------
## Inputs
## ----------------
##
## INPUT: /var/folders/_1/tkj8dh2n73s3c5fp80h0td380000gn/T//RtmpYVQfSd/filec1ea341d06a1/filec1ea2804cf5a.gpkg
## OUTPUT_HTML_FILE: /var/folders/_1/tkj8dh2n73s3c5fp80h0td380000gn/T//RtmpYVQfSd/grocnn.html
##
##
## 0...10...20...30...40...50...60...70...80...90...100 - done.
##
## ----------------
## Results
## ----------------
##
## EXPECTED_MD: 4572.838994129232
## NN_INDEX: 0.5339201195487651
## OBSERVED_MD: 2441.530742422734
## OUTPUT_HTML_FILE: /var/folders/_1/tkj8dh2n73s3c5fp80h0td380000gn/T//RtmpYVQfSd/grocnn.html
## POINT_COUNT: 349
## Z_SCORE: -16.657274867190246