library(mapview)
library(sf)
library(censusxy)
library(dplyr)
library(tidycensus)
library(maptools)Point Pattern Analysis; Drinking Locations in San Antonio, TX
When observing drinking locations in the West side of San Antonio, a Nearest Neighborhoods Index of 1.03 and a Z-score of .29 shows that there is no statistically significant clustering or dispersion within this area. Being that the locations observed are due to randomness with no true differences found among the spacing or clustering of those drinking locations in San Antonio.
From the Kernal Desnsity plots it observes some clustering of drinking locations within the downtown area. When the drinking locations are coupled with the percentage of those living in poverty in San Antonio, only those drinking locations within the downtown area are exhibiting low levels of poverty. As drinking locations move further west the percentage of those living in poverty around those drinking locations increases.
addr2<-read.csv("C:/Users/BTP/Documents/GIS CLASS/Lab 4 HW 5/drinkingSA.csv")
addr2<-addr2[c(4, 6:9, 13:16)]
names(addr2)<-c("name", "street", "city", "st", "zip","employees", "revenue","home based","franchise")
head(addr2) name street city st zip employees
1 Station Cafe 108 King William San Antonio TX 78204 25
2 Splash Bar 905 Nogalitos St San Antonio TX 78204 23
3 Kozy Korner Saloon 507 Ruiz St San Antonio TX 78207 23
4 El Mercado Snack Bar 514 W Commerce St San Antonio TX 78207 19
5 Maggies Tamales & Barbacoa 5721 W Commerce St San Antonio TX 78237 19
6 Pit Stop Bar-B-Q Plus 254 Hobart St San Antonio TX 78237 19
revenue home based franchise
1 $1,756,000 No NA
2 $1,717,000 No NA
3 $1,515,000 No NA
4 $1,234,000 No NA
5 $1,234,000 No NA
6 $1,234,000 No NAresults1<-cxy_geocode(addr2,
street = "street",
city = "city",
state ="st",
zip = "zip",
class="sf",
output = "simple")
results.proj2<-st_transform(results1,
crs = 2278)mapview(results.proj2)mean feature - average of coordinates
mean_feature2<-apply(st_coordinates(results.proj2), MARGIN = 2, FUN = mean)
mean_feature2<-data.frame(place="meanfeature", x=mean_feature2[1], y= mean_feature2[2])
mean_feature2<-st_as_sf(mean_feature2, coords = c("x", "y"), crs= 2278)
mapview(mean_feature2, col.regions="red")+mapview( results1)Central feature - Median of coordinates
median_feature2<-apply(st_coordinates(results.proj2), MARGIN = 2, FUN = median)
median_feature2<-data.frame(place="medianfeature", x=median_feature2[1], y= median_feature2[2])
median_feature2<-st_as_sf(median_feature2, coords = c("x", "y"), crs= 2278)
mapview(median_feature2, col.regions="green")+
mapview(mean_feature2, col.regions="red")+
mapview( results1)Buffer points
drinkingbuff<- st_buffer(results.proj2, dist = 2500)
mapview(drinkingbuff)+mapview(results.proj2, col.regions="green")Convex hull plot
chulldrinking <- st_convex_hull(st_union(results1))
mapview(chulldrinking)+
mapview(results1, 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\BTP\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.batTrying command 'C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat'Success!QGIS version: 3.28.2-FirenzeSaving configuration to 'C:\Users\BTP\AppData\Local/R-qgisprocess/R-qgisprocess/Cache/cache-0.0.0.9000.rds'Metadata of 1170 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 ),]# A tibble: 8 × 27
provider provi…¹ algor…² algor…³ algor…⁴ provi…⁵ can_c…⁶ depre…⁷ group has_k…⁸
<chr> <chr> <chr> <chr> <chr> <chr> <lgl> <lgl> <chr> <lgl>
1 grass7 GRASS grass7… r.li.e… r.li.e… grass7 TRUE FALSE Rast… FALSE
2 grass7 GRASS grass7… r.li.e… r.li.e… grass7 TRUE FALSE Rast… FALSE
3 grass7 GRASS grass7… r.li.p… r.li.p… grass7 TRUE FALSE Rast… FALSE
4 grass7 GRASS grass7… r.li.p… r.li.p… grass7 TRUE FALSE Rast… FALSE
5 native QGIS (… native… linede… Line d… native TRUE FALSE Inte… FALSE
6 qgis QGIS qgis:h… heatma… Heatma… qgis TRUE FALSE Inte… FALSE
7 saga SAGA saga:f… fragme… Fragme… saga TRUE FALSE Rast… FALSE
8 saga SAGA saga:k… kernel… Kernel… saga TRUE FALSE Rast… 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>, …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>
HeatmapRun the algorithm
drinking_dens<-qgis_run_algorithm(algorithm ="qgis:heatmapkerneldensityestimation",
INPUT=results.proj2,
RADIUS = 5000,
PIXEL_SIZE = 100,
KERNEL = 0,
OUTPUT=file.path(getwd(), "drinkingdenst.TIF"),
load_output = TRUE)JSON input ----
{
"inputs": {
"INPUT": "C:\\Users\\BTP\\AppData\\Local\\Temp\\RtmpMBcCSb\\file3a7812076d80\\file3a787b284d20.gpkg",
"RADIUS": 5000,
"PIXEL_SIZE": 100,
"KERNEL": 0,
"OUTPUT_VALUE": 0,
"OUTPUT": "C:/Users/BTP/Documents/GIS CLASS/Lab 4 HW 5/drinkingdenst.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)
resultQ<- qgis_as_raster(drinking_dens)
projection(resultQ)<-crs(results.proj2)
mapview(resultQ)+mapview(results.proj2)Spatial join
A spatial join can combine attributes of one layer with another layer. Here I combine census variables with the Drinking locations 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 ACSDownloading 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.proj2, sa_acs2)
head(spjoin)Simple feature collection with 6 features and 17 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 2100089 ymin: 13694870 xmax: 2130647 ymax: 13705940
Projected CRS: NAD83 / Texas South Central (ftUS)
name street city st zip employees
2 Station Cafe 108 King William San Antonio TX 78204 25
23 Splash Bar 905 Nogalitos St San Antonio TX 78204 23
12 Kozy Korner Saloon 507 Ruiz St San Antonio TX 78207 23
14 El Mercado Snack Bar 514 W Commerce St San Antonio TX 78207 19
15 Maggies Tamales & Barbacoa 5721 W Commerce St San Antonio TX 78237 19
9 Pit Stop Bar-B-Q Plus 254 Hobart St San Antonio TX 78237 19
revenue home based franchise GEOID totpop pblack pwhite phisp
2 $1,756,000 No NA 48029192100 2299 0.7 42.4 0.0
23 $1,717,000 No NA 48029150100 5588 0.8 83.1 2.3
12 $1,515,000 No NA 48029110600 5293 1.1 65.0 6.2
14 $1,234,000 No NA 48029110100 3160 0.5 52.0 0.4
15 $1,234,000 No NA 48029171200 4065 4.6 89.1 3.2
9 $1,234,000 No NA 48029171501 2609 0.0 77.5 0.2
punemp medhhinc ppov geometry
2 3.8 73194 1.2 POINT (2130647 13699644)
23 7.8 40302 14.5 POINT (2124919 13694869)
12 16.1 15250 37.8 POINT (2125077 13705944)
14 2.0 50865 24.7 POINT (2128147 13702613)
15 4.0 28284 22.3 POINT (2104639 13705051)
9 6.6 35194 23.6 POINT (2100089 13698413)mapview(spjoin["ppov"])+mapview(sa_trol)library(tmap)
library(tmaptools)
map5 <- tm_shape(sa_acs2)+
tm_polygons()+
tm_shape(spjoin)+
tm_dots("ppov", title="% in Poverty",
palette="Reds",
style="pretty",
n=5,
size=0.2)+
tm_format("World",
main.title="Poverty at Drinking Locations, San Anonio's Westside",
main.title.position=c('center','top'),
main.title.size=1.5,
title="Author: B.A. Flores, M.S. \nSource: ACS 2019 \nDetails:Poverty estimates around Westside San Antonio Drinking Locations",
legend.title.size=1.7,
legend.outside=T,
legend.text.size=1.2)+
tm_scale_bar(position = c("left","bottom"))+
tm_compass()
map5
Nearest Neighbor analysis
library(spatstat)Loading required package: spatstat.dataLoading required package: spatstat.geomspatstat.geom 3.0-6
Attaching package: 'spatstat.geom'The following objects are masked from 'package:raster':
area, rotate, shiftLoading required package: spatstat.randomspatstat.random 3.1-3Loading required package: spatstat.exploreLoading required package: nlme
Attaching package: 'nlme'The following object is masked from 'package:raster':
getDataThe following object is masked from 'package:dplyr':
collapsespatstat.explore 3.0-6Loading required package: spatstat.modelLoading required package: rpartspatstat.model 3.2-1Loading required package: spatstat.linnetspatstat.linnet 3.0-4
spatstat 3.0-3
For an introduction to spatstat, type 'beginner' drinking.pp<-as.ppp(as(results.proj2, "Spatial"))
plot(nearest.neighbour(drinking.pp))
algs[grepl(pattern = "nearest", x = algs$algorithm ),]# A tibble: 9 × 27
provider provi…¹ algor…² algor…³ algor…⁴ provi…⁵ can_c…⁶ depre…⁷ group has_k…⁸
<chr> <chr> <chr> <chr> <chr> <chr> <lgl> <lgl> <chr> <lgl>
1 gdal GDAL gdal:g… gridin… Grid (… gdal TRUE FALSE Rast… FALSE
2 gdal GDAL gdal:g… gridne… Grid (… gdal TRUE FALSE Rast… FALSE
3 native QGIS (… native… anglet… Align … native TRUE FALSE Cart… FALSE
4 native QGIS (… native… joinby… Join a… native TRUE FALSE Vect… FALSE
5 native QGIS (… native… neares… Neares… native TRUE FALSE Vect… FALSE
6 qgis QGIS qgis:d… distan… Distan… qgis TRUE FALSE Vect… FALSE
7 qgis QGIS qgis:d… distan… Distan… qgis TRUE FALSE Vect… FALSE
8 qgis QGIS qgis:k… kneare… Concav… qgis TRUE FALSE Vect… FALSE
9 saga SAGA saga:n… neares… Neares… saga TRUE FALSE Rast… 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>, …qgis_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 (optional)
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-scoredrinking_nn<-qgis_run_algorithm(alg="native:nearestneighbouranalysis",
INPUT=results.proj2,
OUTPUT_HTML_FILE=file.path(tempdir(), "drinkingnn.html"),
load_output = TRUE)Ignoring unknown input 'load_output'JSON input ----
{
"inputs": {
"INPUT": "C:\\Users\\BTP\\AppData\\Local\\Temp\\RtmpMBcCSb\\file3a7812076d80\\file3a78375e389f.gpkg",
"OUTPUT_HTML_FILE": "C:\\Users\\BTP\\AppData\\Local\\Temp\\RtmpMBcCSb/drinkingnn.html"
}
}
Running cmd.exe /c call \
"C:/Program Files/QGIS 3.28.2/bin/qgis_process-qgis.bat" --json run \
"native:nearestneighbouranalysis" -drinking_nn<Result of `qgis_run_algorithm("native:nearestneighbouranalysis", ...)`>
List of 6
$ EXPECTED_MD : num 2045
$ NN_INDEX : num 1.03
$ OBSERVED_MD : num 2113
$ OUTPUT_HTML_FILE: 'qgis_outputHtml' chr "C:\\Users\\BTP\\AppData\\Local\\Temp\\RtmpMBcCSb/drinkingnn.html"
$ POINT_COUNT : num 21
$ Z_SCORE : num 0.29