*These are updates. The 2021 version from Corey Sparks, Ph.D. can be found here. The 2020 version is here. The 2018 version is here.

This example shows how to use R and QGIS from within R to perform a series of common point pattern analysis techniques.

library(mapview)
library(sf)
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(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, union
addr<-read.csv(url("https://raw.githubusercontent.com/coreysparks/DEM7093/main/data/west_side_groceries.csv"))
addr<-addr[c(6, 12:14)]
names(addr)<-c("street", "city", "st", "zip")
head(addr)
results<-cxy_geocode(addr,
                     street = "street",
                     city = "city",
                     state ="st",
                     zip = "zip",
                     class="sf",
                     output = "simple")
## 14 rows removed to create an sf object. These were addresses that the geocoder could not match.
results.proj<-st_transform(results,
                           crs = 2278)

OR just use the lat / long information in the data!

addr<-read.csv(url("https://raw.githubusercontent.com/coreysparks/DEM7093/main/data/west_side_groceries.csv"))
results <- st_as_sf(addr, coords=c("Longitude", "Latitude"), crs=4269,agr="constant")
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

wicbuff<- st_buffer(results.proj, dist = 2500)
mapview(wicbuff)+mapview(results.proj, col.regions="green")

Convex hull plot

chull <- st_convex_hull(st_union(results))
mapview(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 '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()
## 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.

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 ),]
qgis_show_help("qgis:heatmapkerneldensityestimation")
## Heatmap (Kernel Density Estimation) (qgis:heatmapkerneldensityestimation)
## 
## ----------------
## Description
## ----------------
## Creates a density (heatmap) raster of an input point vector layer using kernel density estimation. Heatmaps allow easy identification of hotspots and clustering of points.
## The density is calculated based on the number of points in a location, with larger numbers of clustered points resulting in larger values.
## 
## 
## ----------------
## Arguments
## ----------------
## 
## INPUT: Point layer
##  Argument type:  source
##  Acceptable values:
##      - Path to a vector layer
## RADIUS: Radius
##  Default value:  100
##  Argument type:  distance
##  Acceptable values:
##      - A numeric value
## RADIUS_FIELD: Radius from field (optional)
##  Argument type:  field
##  Acceptable values:
##      - The name of an existing field
##      - ; delimited list of existing field names
## PIXEL_SIZE: Output raster size
##  Default value:  0.1
##  Argument type:  number
##  Acceptable values:
##      - A numeric value
## WEIGHT_FIELD: Weight from field (optional)
##  Argument type:  field
##  Acceptable values:
##      - The name of an existing field
##      - ; delimited list of existing field names
## KERNEL: Kernel shape
##  Default value:  0
##  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) (optional)
##  Default value:  0
##  Argument type:  number
##  Acceptable values:
##      - A numeric value
## OUTPUT_VALUE: Output value scaling
##  Default value:  0
##  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>
##  Heatmap

Run the algorithm

wic_dens<-qgis_run_algorithm(algorithm ="qgis:heatmapkerneldensityestimation",
         INPUT=results.proj,
         RADIUS = 5000,
         PIXEL_SIZE = 100,
         KERNEL = 0,
         OUTPUT=file.path(getwd(), "wicdenst.TIF"),
         load_output = TRUE)
## JSON input ----
## {
##   "inputs": {
##     "INPUT": "C:\\Users\\xee291\\AppData\\Local\\Temp\\RtmpkfEJIj\\file1a6066515bc9\\file1a605e641908.gpkg",
##     "RADIUS": 5000,
##     "PIXEL_SIZE": 100,
##     "KERNEL": 0,
##     "OUTPUT_VALUE": 0,
##     "OUTPUT": "C:/Users/xee291/OneDrive - University of Texas at San Antonio/Documents/UTSA 2022-2023/GIS Class/Week 7/wicdenst.TIF"
##   }
## }
## 
## Running cmd.exe /c call \
##   "C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat" --json run \
##   "qgis:heatmapkerneldensityestimation" -
library(raster)
library(RColorBrewer)
result<- qgis_as_raster(wic_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 WIC 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)
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$nwic<- lengths(st_intersects(sa_acs2, results.proj))
mapview(sa_acs2, zcol="nwic")+
  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 ),]
qgis_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)
##  Default value:  0
##  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 polygons
wic_von<-qgis_run_algorithm(alg="qgis:voronoipolygons",
         INPUT=results.proj,
         OUTPUT=file.path(tempdir(), "wicvon.shp"),
         load_output = TRUE)
## JSON input ----
## {
##   "inputs": {
##     "INPUT": "C:\\Users\\xee291\\AppData\\Local\\Temp\\RtmpkfEJIj\\file1a6066515bc9\\file1a603b25fce.gpkg",
##     "OUTPUT": "C:\\Users\\xee291\\AppData\\Local\\Temp\\RtmpkfEJIj/wicvon.shp"
##   }
## }
## 
## Running cmd.exe /c call \
##   "C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat" --json run \
##   "qgis:voronoipolygons" -
wic_von<-sf::read_sf(qgis_output(wic_von, "OUTPUT"))
mapview(wic_von, alpha=.75)+
  mapview(results.proj, col.regions="green")

Nearest Neighbor analysis

  • NOTE: I have not gotten these chunks to run yet. See Dr. Sparks’ here for what the output should look like.
library(spatstat)
wic.pp<-as.ppp(as(results.proj, "Spatial"))
plot(nearest.neighbour(wic.pp))
algs[grepl(pattern = "nearest", x = algs$algorithm ),]
qgis_show_help("native:nearestneighbouranalysis")
wic_nn<-qgis_run_algorithm(alg="native:nearestneighbouranalysis",
         INPUT=results.proj,
        OUTPUT_HTML_FILE=file.path(tempdir(), "wicnn.html"),
         load_output = TRUE)