Density estimation
mean feature

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.6.1, GDAL 2.2.3, PROJ 4.9.3

library(censusxy)

addr<-read.csv(url("https://raw.githubusercontent.com/coreysparks/data/master/wic_west_side.csv"))
addr<-addr[c(6, 12:14)]
names(addr)<-c("street", "city", "st", "zip")
head(addr)

ABCDEFGHIJ0123456789

	street <fctr>	city <fctr>	st <fctr>	zip <int>
1	1039 W Hildebrand Ave	San Antonio	TX	78201
2	1102 Delgado St	San Antonio	TX	78207
3	1115 W Martin St	San Antonio	TX	78207
4	1115 W Martin St	San Antonio	TX	78207
5	1115 W Martin St	San Antonio	TX	78207
6	1115 W Martin St	San Antonio	TX	78207

results<-cxy_geocode(addr, address = street, city = city, state =st, zip = zip, output = "sf")

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 21 rows [5, 6, 8,
## 12, 17, 20, 26, 32, 33, 40, 46, 48, 50, 51, 52, 54, 55, 57, 58, 60, ...].

results.proj<-st_transform(results, crs = 2278)

mapview(results)

Density estimation

library(RQGIS3)

## Loading required package: reticulate

set_env()

## Trying to find QGIS in C:/OSGEO4~1

## $root
## [1] "C:/OSGeo4W64"
## 
## $qgis_prefix_path
## [1] "C:/OSGeo4W64/apps/qgis"
## 
## $python_plugins
## [1] "C:/OSGeo4W64/apps/qgis/python/plugins"
## 
## $platform
## [1] "Windows"

open_app() #takes a while

mean feature

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)

Buffer points

wicbuff<- st_buffer(results.proj, dist = 2500)
mapview(wicbuff)

Convex hull plot

find_algorithms(search_term = "minimum")

## [1] "Minimum enclosing circles---------------------------->native:minimumenclosingcircle"    
## [2] "Oriented minimum bounding box------------------------>native:orientedminimumboundingbox"
## [3] "Minimum bounding geometry---------------------------->qgis:minimumboundinggeometry"     
## [4] "Minimum distance analysis---------------------------->saga:minimumdistanceanalysis"

get_usage(alg="qgis:minimumboundinggeometry")

## Minimum bounding geometry (qgis:minimumboundinggeometry)
## 
## This algorithm creates geometries which enclose the features from an input layer.
## Numerous enclosing geometry types are supported
## including bounding boxes (envelopes)
## oriented rectangles
## circles and convex hulls.
## Optionally
## the features can be grouped by a field. If set
## this causes the output layer to contain one feature per grouped value with a minimal geometry covering just the features with matching values.
## 
## 
## 
## ----------------
## Input parameters
## ----------------
## 
## INPUT: Input layer
## 
##  Parameter type: QgsProcessingParameterFeatureSource
## 
##  Accepted data types:
##      - str: layer ID
##      - str: layer name
##      - str: layer source
##      - QgsProcessingFeatureSourceDefinition
##      - QgsProperty
##      - QgsVectorLayer
## 
## FIELD: Field (optional
## set if features should be grouped by class)
## 
##  Parameter type: QgsProcessingParameterField
## 
##  Accepted data types:
##      - str
##      - QgsProperty
## 
## TYPE: Geometry type
## 
##  Parameter type: QgsProcessingParameterEnum
## 
##  Available values:
##      - 0: Envelope (Bounding Box)
##      - 1: Minimum Oriented Rectangle
##      - 2: Minimum Enclosing Circle
##      - 3: Convex Hull
## 
##  Accepted data types:
##      - int
##      - str: as string representation of int
## e.g. 1
##      - QgsProperty
## 
## OUTPUT: Bounding geometry
## 
##  Parameter type: QgsProcessingParameterFeatureSink
## 
##  Accepted data types:
##      - str: destination vector file
## e.g. d:/test.shp
##      - str: memory: to store result in temporary memory layer
##      - str: using vector provider ID prefix and destination URI
## e.g. postgres:... to store result in PostGIS table
##      - QgsProcessingOutputLayerDefinition
##      - QgsProperty
## 
## ----------------
## Outputs
## ----------------
## 
## OUTPUT:  <QgsProcessingOutputVectorLayer>
##  Bounding geometry

wic_hull<-run_qgis(alg="qgis:minimumboundinggeometry",
         INPUT=results.proj,
         TYPE=3,
         OUTPUT=file.path(tempdir(), "wichull.shp"),
         load_output = TRUE)

## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver

## $OUTPUT
## [1] "C:/Users/ozd504/AppData/Local/Temp/RtmpyaTzTk/wichull.shp"

mapview(wic_hull)+mapview(results.proj)

kernel density - You need projected data for this to work right

find_algorithms(search_term = "density")

## [1] "r.li.edgedensity------------------------------------->grass7:r.li.edgedensity"                             
## [2] "r.li.edgedensity.ascii------------------------------->grass7:r.li.edgedensity.ascii"                       
## [3] "r.li.patchdensity------------------------------------>grass7:r.li.patchdensity"                            
## [4] "r.li.patchdensity.ascii------------------------------>grass7:r.li.patchdensity.ascii"                      
## [5] "Heatmap (Kernel Density Estimation)------------------>qgis:heatmapkerneldensityestimation"                 
## [6] "Fragmentation classes from density and connectivity--->saga:fragmentationclassesfromdensityandconnectivity"
## [7] "Kernel density estimation---------------------------->saga:kerneldensityestimation"

get_usage(alg="qgis:heatmapkerneldensityestimation")

## Heatmap (Kernel Density Estimation) (qgis:heatmapkerneldensityestimation)
## 
## 
## ----------------
## Input parameters
## ----------------
## 
## INPUT: Point layer
## 
##  Parameter type: QgsProcessingParameterFeatureSource
## 
##  Accepted data types:
##      - str: layer ID
##      - str: layer name
##      - str: layer source
##      - QgsProcessingFeatureSourceDefinition
##      - QgsProperty
##      - QgsVectorLayer
## 
## RADIUS: Radius
## 
##  Parameter type: QgsProcessingParameterDistance
## 
##  Accepted data types:
##      - int
##      - float
##      - QgsProperty
## 
## RADIUS_FIELD: Radius from field
## 
##  Parameter type: QgsProcessingParameterField
## 
##  Accepted data types:
##      - str
##      - QgsProperty
## 
## PIXEL_SIZE: Output raster size
## 
##  Parameter type: ParameterHeatmapPixelSize
## 
##  Accepted data types:
##      - int
##      - float
##      - QgsProperty
## 
## WEIGHT_FIELD: Weight from field
## 
##  Parameter type: QgsProcessingParameterField
## 
##  Accepted data types:
##      - str
##      - QgsProperty
## 
## KERNEL: Kernel shape
## 
##  Parameter type: QgsProcessingParameterEnum
## 
##  Available values:
##      - 0: Quartic
##      - 1: Triangular
##      - 2: Uniform
##      - 3: Triweight
##      - 4: Epanechnikov
## 
##  Accepted data types:
##      - int
##      - str: as string representation of int
## e.g. 1
##      - QgsProperty
## 
## DECAY: Decay ratio (Triangular kernels only)
## 
##  Parameter type: QgsProcessingParameterNumber
## 
##  Accepted data types:
##      - int
##      - float
##      - QgsProperty
## 
## OUTPUT_VALUE: Output value scaling
## 
##  Parameter type: QgsProcessingParameterEnum
## 
##  Available values:
##      - 0: Raw
##      - 1: Scaled
## 
##  Accepted data types:
##      - int
##      - str: as string representation of int
## e.g. 1
##      - QgsProperty
## 
## OUTPUT: Heatmap
## 
##  Parameter type: QgsProcessingParameterRasterDestination
## 
##  Accepted data types:
##      - str
##      - QgsProperty
##      - QgsProcessingOutputLayerDefinition
## 
## ----------------
## Outputs
## ----------------
## 
## OUTPUT:  <QgsProcessingOutputRasterLayer>
##  Heatmap

wic_dens<-run_qgis(alg="qgis:heatmapkerneldensityestimation",
         INPUT=results.proj,
         RADIUS = 5000,
         PIXEL_SIZE = 100,
         KERNEL = 0,
         OUTPUT=file.path(tempdir(), "wichull.TIF"),
         load_output = TRUE)

## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver

## $OUTPUT
## [1] "C:/Users/ozd504/AppData/Local/Temp/RtmpyaTzTk/wichull.TIF"

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

library(raster)

## Loading required package: sp

library(RColorBrewer)
mapview(wic_dens, alpha=.75)+mapview(results.proj)

## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files

Spatial join

library(tidycensus)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#load census tract data
sa_acs<-get_acs(geography = "tract", state="TX", county = c("Bexar"), year = 2018,
                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 2014-2018 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, pfemhh=DP02_0008PE,
         phsormore=DP02_0066PE,punemp=DP03_0009PE, medhhinc=DP03_0062E,
         ppov=DP03_0119PE)%>%
  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)

ABCDEFGHIJ0123456789

	street <fctr>	city <fctr>	st <fctr>	zip <int>	cxy_status <chr>	cxy_quality <chr>
1	1039 W Hildebrand Ave	San Antonio	TX	78201	Match	Non_Exact
2	1102 Delgado St	San Antonio	TX	78207	Match	Exact
3	1115 W Martin St	San Antonio	TX	78207	Match	Non_Exact
4	1115 W Martin St	San Antonio	TX	78207	Match	Non_Exact
5	1115 W Martin St	San Antonio	TX	78207	Match	Non_Exact
6	1115 W Martin St	San Antonio	TX	78207	Match	Non_Exact

mapview(spjoin["punemp"])+mapview(sa_trol)

Thiessen/Voronoi Polygons

find_algorithms(search_term = "voronoi")

## [1] "v.voronoi-------------------------------------------->grass7:v.voronoi"    
## [2] "Voronoi polygons------------------------------------->qgis:voronoipolygons"

get_usage(alg="qgis:voronoipolygons")

## Voronoi polygons (qgis:voronoipolygons)
## 
## This algorithm takes a points layer and generates a polygon layer containing the voronoi polygons corresponding to those input points.
## 
## 
## 
## ----------------
## Input parameters
## ----------------
## 
## INPUT: Input layer
## 
##  Parameter type: QgsProcessingParameterFeatureSource
## 
##  Accepted data types:
##      - str: layer ID
##      - str: layer name
##      - str: layer source
##      - QgsProcessingFeatureSourceDefinition
##      - QgsProperty
##      - QgsVectorLayer
## 
## BUFFER: Buffer region (% of extent)
## 
##  Parameter type: QgsProcessingParameterNumber
## 
##  Accepted data types:
##      - int
##      - float
##      - QgsProperty
## 
## OUTPUT: Voronoi polygons
## 
##  Parameter type: QgsProcessingParameterFeatureSink
## 
##  Accepted data types:
##      - str: destination vector file
## e.g. d:/test.shp
##      - str: memory: to store result in temporary memory layer
##      - str: using vector provider ID prefix and destination URI
## e.g. postgres:... to store result in PostGIS table
##      - QgsProcessingOutputLayerDefinition
##      - QgsProperty
## 
## ----------------
## Outputs
## ----------------
## 
## OUTPUT:  <QgsProcessingOutputVectorLayer>
##  Voronoi polygons

wic_von<-run_qgis(alg="qgis:voronoipolygons",
         INPUT=results.proj,
         OUTPUT=file.path(tempdir(), "wicvon.shp"),
         load_output = TRUE)

## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver

## $OUTPUT
## [1] "C:/Users/ozd504/AppData/Local/Temp/RtmpyaTzTk/wicvon.shp"

mapview(wic_von, alpha=.75)+mapview(results.proj)

Nearest Neighbor analysis

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## The following object is masked from 'package:raster':
## 
##     getData

## Loading required package: rpart

## 
## spatstat 1.61-0       (nickname: 'Puppy zoomies') 
## For an introduction to spatstat, type 'beginner'

## 
## Note: spatstat version 1.61-0 is out of date by more than 5 months; we recommend upgrading to the latest version.

## 
## Attaching package: 'spatstat'

## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift

wic.pp<-as.ppp(results.proj)

## Warning in as.ppp.sf(results.proj): only first attribute column is used for
## marks

plot(nearest.neighbour(wic.pp))

find_algorithms(search_term = "analysis")

##  [1] "Nearest neighbour analysis--------------------------->qgis:nearestneighbouranalysis"                                           
##  [2] "Basic terrain analysis------------------------------->saga:basicterrainanalysis"                                               
##  [3] "Cluster analysis------------------------------------->saga:clusteranalysis"                                                    
##  [4] "Cluster analysis (shapes)---------------------------->saga:clusteranalysisshapes"                                              
##  [5] "Connectivity analysis-------------------------------->saga:connectivityanalysis"                                               
##  [6] "Fire risk analysis----------------------------------->saga:fireriskanalysis"                                                   
##  [7] "Minimum distance analysis---------------------------->saga:minimumdistanceanalysis"                                            
##  [8] "Multiple linear regression analysis------------------>saga:multiplelinearregressionanalysis"                                   
##  [9] "Multiple linear regression analysis (shapes)--------->saga:multiplelinearregressionanalysisshapes"                             
## [10] "Multiple regression analysis (grid and predictor grids)--->saga:multipleregressionanalysisgridandpredictorgrids"               
## [11] "Multiple regression analysis (points and predictor grids)--->saga:multipleregressionanalysispointsandpredictorgrids"           
## [12] "Multiple regression analysis (points/raster)--------->saga:multipleregressionanalysispointsraster"                             
## [13] "Multiple regression analysis (raster/raster)--------->saga:multipleregressionanalysisrasterraster"                             
## [14] "Pattern analysis------------------------------------->saga:patternanalysis"                                                    
## [15] "Principle components analysis------------------------>saga:principlecomponentsanalysis"                                        
## [16] "Regression analysis---------------------------------->saga:regressionanalysis"                                                 
## [17] "Regression analysis (points and predictor grid)------>saga:regressionanalysispointsandpredictorgrid"                           
## [18] "Residual analysis------------------------------------>saga:residualanalysis"                                                   
## [19] "Spatial point pattern analysis----------------------->saga:spatialpointpatternanalysis"                                        
## [20] "Zonal multiple regression analysis (points and predictor grids)--->saga:zonalmultipleregressionanalysispointsandpredictorgrids"

get_usage(alg="qgis:nearestneighbouranalysis")

## Nearest neighbour analysis (qgis:nearestneighbouranalysis)
## 
## This algorithm performs nearest neighbor analysis for a point layer.
## Output is generated as an html file with the computed statistical values.
## 
## 
## 
## ----------------
## Input parameters
## ----------------
## 
## INPUT: Input layer
## 
##  Parameter type: QgsProcessingParameterFeatureSource
## 
##  Accepted data types:
##      - str: layer ID
##      - str: layer name
##      - str: layer source
##      - QgsProcessingFeatureSourceDefinition
##      - QgsProperty
##      - QgsVectorLayer
## 
## OUTPUT_HTML_FILE: Nearest neighbour
## 
##  Parameter type: QgsProcessingParameterFileDestination
## 
##  Accepted data types:
##      - str
##      - QgsProperty
## 
## ----------------
## Outputs
## ----------------
## 
## OUTPUT_HTML_FILE:  <QgsProcessingOutputHtml>
##  Nearest neighbour
## 
## OBSERVED_MD:  <QgsProcessingOutputNumber>
##  Observed mean distance
## 
## EXPECTED_MD:  <QgsProcessingOutputNumber>
##  Expected mean distance
## 
## NN_INDEX:  <QgsProcessingOutputNumber>
##  Nearest neighbour index
## 
## POINT_COUNT:  <QgsProcessingOutputNumber>
##  Number of points
## 
## Z_SCORE:  <QgsProcessingOutputNumber>
##  Z-Score

wic_nn<-run_qgis(alg="qgis:nearestneighbouranalysis",
         INPUT=results.proj,
        OUTPUT_HTML_FILE=file.path(tempdir(), "wicnn.html"))

## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver

## $EXPECTED_MD
## [1] 1700.572
## 
## $NN_INDEX
## [1] 0.5277587
## 
## $OBSERVED_MD
## [1] 897.4918
## 
## $OUTPUT_HTML_FILE
## [1] "C:/Users/ozd504/AppData/Local/Temp/RtmpyaTzTk/wicnn.html"
## 
## $POINT_COUNT
## [1] 62
## 
## $Z_SCORE
## [1] -7.11362

DEM 5093/7093 Point Pattern Analysis

Corey Sparks, Ph.D.

2/27/2020

Density estimation

mean feature

Buffer points

Convex hull plot

kernel density - You need projected data for this to work right

Spatial join

Thiessen/Voronoi Polygons

Nearest Neighbor analysis