Dem5093_HW4
Sarah Sharmin
2/27/2020
Objectives
For the data on west side grocery stores, show:
1.2500 foot buffers
2. Convex hull of the points
3. Spatial join of the points to the ACS layer
4. Nearest Neighbor analysis of the grocery store layer. Report the numeric summary of the analysis.
Loading Library
library(mapview)
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3library(censusxy)Loading west side grocery store data
addr<-read.csv("C:/Users/sarah/Google Drive/MSc Demography/Spring 2020/GIS 5093/gis_class/Homework3/dem5093_Gstore.csv")
addr<-addr[c(6, 12:14)]
names(addr)<-c("street", "city", "st", "zip")
head(addr)## street city st zip
## 1 1545 S San Marcos San Antonio TX 78207
## 2 4522 Fredericksburg Rd San Antonio TX 78201
## 3 8234 Agora Pkwy San Antonio TX 78201
## 4 2902 Guadalupe St San Antonio TX 78207
## 5 1200 El Paso St San Antonio TX 78207
## 6 516 S Flores St San Antonio TX 78204results<-cxy_geocode(addr, address = street, city = city, state =st, zip = zip, output = "sf")## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 24 rows [1, 5, 7,
## 9, 14, 21, 22, 26, 27, 28, 29, 31, 34, 43, 45, 46, 50, 51, 53, 54, ...].results.proj<-st_transform(results, crs = 2278)A point Map view of the west side grocery store
mapview(results)Density estimation
library(RQGIS3)## Loading required package: reticulateset_env()## Trying to find QGIS in C:/## $root
## [1] "C:/Program Files/QGIS 3.4"
##
## $qgis_prefix_path
## [1] "C:/Program Files/QGIS 3.4/apps/qgis-ltr"
##
## $python_plugins
## [1] "C:/Program Files/QGIS 3.4/apps/qgis-ltr/python/plugins"
##
## $platform
## [1] "Windows"open_app() #takes a whilemean feature
To see where the mean point is located.
mean_feature<-apply(st_coordinates(results.proj), MARGIN = 2, FUN = median)
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
To see coverage area of each grocery store. Here we used 2500 feet point buffer
gsbuff<- st_buffer(results.proj, dist = 2500)
mapview(gsbuff)Convex hull plot
To see the overall spatial extent of all the grocery stores together.
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 geometrygs_hull<-run_qgis(alg="qgis:minimumboundinggeometry",
INPUT=results.proj,
TYPE=3,
OUTPUT=file.path(tempdir(), "gshull.shp"),
load_output = TRUE)## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver## $OUTPUT
## [1] "C:/Users/sarah/AppData/Local/Temp/RtmpqGbRJm/gshull.shp"mapview(gs_hull)+mapview(results.proj)kernel density
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>
## Heatmapgs_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/sarah/AppData/Local/Temp/RtmpqGbRJm/wichull.TIF"## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared fileslibrary(raster)## Loading required package: splibrary(RColorBrewer)
mapview(gs_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 filesSpatial join
joining coordinates with ACS 2018 data
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)## Simple feature collection with 6 features and 15 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 2116747 ymin: 13700040 xmax: 2128732 ymax: 13709460
## epsg (SRID): 2278
## proj4string: +proj=lcc +lat_1=30.28333333333333 +lat_2=28.38333333333333 +lat_0=27.83333333333333 +lon_0=-99 +x_0=600000 +y_0=3999999.9998984 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
## street city st zip cxy_status cxy_quality
## 1 2902 Guadalupe St San Antonio TX 78207 Match Exact
## 2 516 S Flores St San Antonio TX 78204 Match Non_Exact
## 3 2614 W Commerce St San Antonio TX 78207 Match Non_Exact
## 4 1557 Sw 19th St San Antonio TX 78207 Match Non_Exact
## 5 1424 Guadalupe St Ste 2 San Antonio TX 78207 Match Exact
## 6 326 N Zarzamora St San Antonio TX 78207 Match Non_Exact
## cxy_match GEOID totpop pblack pwhite
## 1 2902 GUADALUPE ST, SAN ANTONIO, TX, 78207 48029170300 6350 0.1 87.2
## 2 516 S FLORES ST, SAN ANTONIO, TX, 78204 48029192100 2064 1.1 39.6
## 3 2614 W COMMERCE ST, SAN ANTONIO, TX, 78207 48029170200 5507 0.4 92.6
## 4 1557 NW 19TH ST, SAN ANTONIO, TX, 78207 48029170401 4358 0.0 76.8
## 5 1424 GUADALUPE ST, SAN ANTONIO, TX, 78207 48029170200 5507 0.4 92.6
## 6 326 S ZARZAMORA ST, SAN ANTONIO, TX, 78207 48029170200 5507 0.4 92.6
## phisp punemp medhhinc ppov geometry
## 1 0.0 10.2 25986 32.0 POINT (2117258 13700560)
## 2 0.0 3.7 68333 1.5 POINT (2128732 13700439)
## 3 1.9 12.2 25691 36.5 POINT (2120109 13703408)
## 4 0.0 11.1 22529 43.8 POINT (2116747 13709457)
## 5 1.9 12.2 25691 36.5 POINT (2122497 13700037)
## 6 1.9 12.2 25691 36.5 POINT (2118991 13702620)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 polygonswic_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/sarah/AppData/Local/Temp/RtmpqGbRJm/wicvon.shp"mapview(wic_von, alpha=.75)+mapview(results.proj)Nearest Neighbor analysis
To see the pattern is clustered or disperesed
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.63-0 (nickname: 'Space camouflage')
## For an introduction to spatstat, type 'beginner'##
## Attaching package: 'spatstat'## The following objects are masked from 'package:raster':
##
## area, rotate, shiftgs.pp<-as.ppp(results.proj)## Warning in as.ppp.sf(results.proj): only first attribute column is used for
## marksplot(nearest.neighbour(gs.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-Scoregs_nn<-run_qgis(alg="qgis:nearestneighbouranalysis",
INPUT=results.proj,
OUTPUT_HTML_FILE=file.path(tempdir(), "gsnn.html"))## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver## $EXPECTED_MD
## [1] 2035.029
##
## $NN_INDEX
## [1] 0.7066506
##
## $OBSERVED_MD
## [1] 1438.054
##
## $OUTPUT_HTML_FILE
## [1] "C:/Users/sarah/AppData/Local/Temp/RtmpqGbRJm/gsnn.html"
##
## $POINT_COUNT
## [1] 38
##
## $Z_SCORE
## [1] -3.459457As the neighborhood index is less than 1, the pattern is clustered.