The goal of this task is to illustrate how to use R and QGIS from within R to perform some of common point pattern analysis techniques.
library(mapview)
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.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
Create Point layer of Bars from A to Z database
addr<-read.csv("C:\\Users\\anami\\OneDrive\\Desktop\\GIS_CLASSFINAL\\drinking_Places.csv")
addr <- addr%>%
dplyr::select(names(addr)[c(4,6,14,23:24,41,42)])
results <- st_as_sf(addr, coords=c("Longitude", "Latitude"), crs=4269,agr="constant")
results.proj<-st_transform(results,
crs = 2278)
mv1 <- mapview(results.proj)
mapshot(mv1, file = paste0(getwd(), "/map1.png"))
## PhantomJS not found. You can install it with webshot::install_phantomjs(). If it is installed, please make sure the phantomjs executable can be found via the PATH variable.
mv1
library(tmap)
library(tmaptools)
library(OpenStreetMap)
## Warning: package 'OpenStreetMap' was built under R version 4.1.3
bg<- read_osm(results.proj, ext=1.1)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots()
tmap_mode("plot")
## tmap mode set to plotting
library(ggplot2)
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
Mean features-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)
tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj, size = .2)+
tm_dots(col = "green", size=.15)+
tm_shape(mean_feature)+
tm_dots(col = "red", size = .4)
Convex hull plot
chull <- st_convex_hull(st_union(results))
tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col = "green", size = .1)+
tm_shape(chull)+
tm_polygons(col = "grey", alpha = .5)
Spatial join
library(tidycensus)
## Warning: package 'tidycensus' was built under R version 4.1.3
sa_acs<-get_acs(geography = "tract",
state="TX",
county = "Bexar",
year = 2019,
variables=c( "DP05_0001E", "DP03_0062E",
"DP03_0119PE","DP02_0066PE",
"DP03_0062E","DP03_0119PE",
"DP05_0073PE","DP05_0066PE") ,
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
sa_acs2<-sa_acs%>%
mutate(totpop= DP05_0001E,
pblack=DP05_0073PE,
phisp=DP05_0066PE,
phsormore=DP02_0066PE,
medhhinc=DP03_0062E,
ppov=DP03_0119PE)%>%
dplyr::select(GEOID, totpop, pblack, phisp, 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)
tmap_mode("view")
## tmap mode set to interactive viewing
Count points in polygons
sa_acs2$nbars<- lengths(st_intersects(sa_acs2, results.proj))
sa_acs2$bars_pc <- 1000*(sa_acs2$nbars/sa_acs2$totpop)
tmap_mode("plot")
## tmap mode set to plotting
m<-tm_basemap("OpenStreetMap.Mapnik")+
tm_shape(sa_acs2)+
tm_polygons("bars_pc")+
tm_shape(spjoin, is.master = T)+
tm_dots( col="blue", size = .4)
m
Nearest neighbor analysis
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.3-1
## Loading required package: spatstat.core
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rpart
## spatstat.core 2.3-2
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-1
##
## spatstat 2.3-0 (nickname: 'That's not important right now')
## For an introduction to spatstat, type 'beginner'
bars.pp<-as.ppp(as(results.proj, "Spatial"))
plot(nearest.neighbour(bars.pp))
options(qgisprocess.path = "C://Program Files//QGIS 3.16.16//bin//qgis_process-qgis-ltr.bat")
library(qgisprocess)
## Using 'qgis_process' at 'C://Program Files//QGIS 3.16.16//bin//qgis_process-qgis-ltr.bat'.
## QGIS version: 3.16.16-Hannover
## Configuration loaded from 'C:\Users\anami\AppData\Local/R-qgisprocess/R-qgisprocess/Cache/cache-0.0.0.9000.rds'
## Run `qgis_configure()` for details.
## Using JSON for output serialization
qgis_configure()
## Trying getOption('qgisprocess.path'): 'C://Program Files//QGIS 3.16.16//bin//qgis_process-qgis-ltr.bat'
## Success!
## QGIS version: 3.16.16-Hannover
## Saving configuration to 'C:\Users\anami\AppData\Local/R-qgisprocess/R-qgisprocess/Cache/cache-0.0.0.9000.rds'
## Metadata of 1134 algorithms queried and stored in cache.
## Run `qgis_algorithms()` to see them.
## - Using JSON for output serialization
algs<-qgis_algorithms()
algs[grepl(pattern = "nearest", x = algs$algorithm ),]
## Registered S3 method overwritten by 'cli':
## method from
## print.boxx spatstat.geom
## # A tibble: 9 x 27
## provider provider_title algorithm algorithm_id algorithm_title provider_id
## <chr> <chr> <chr> <chr> <chr> <chr>
## 1 gdal GDAL gdal:grid~ gridinverse~ Grid (IDW with~ gdal
## 2 gdal GDAL gdal:grid~ gridnearest~ Grid (Nearest ~ gdal
## 3 native QGIS (native c++) native:an~ angletonear~ Align points t~ native
## 4 native QGIS (native c++) native:jo~ joinbyneare~ Join attribute~ native
## 5 native QGIS (native c++) native:ne~ nearestneig~ Nearest neighb~ native
## 6 qgis QGIS qgis:dist~ distanceton~ Distance to ne~ qgis
## 7 qgis QGIS qgis:dist~ distanceton~ Distance to ne~ qgis
## 8 qgis QGIS qgis:knea~ knearestcon~ Concave hull (~ qgis
## 9 saga SAGA saga:near~ nearestneig~ Nearest Neighb~ saga
## # ... with 21 more variables: can_cancel <lgl>, deprecated <lgl>, group <chr>,
## # has_known_issues <lgl>, 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>, ...
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
## 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-score
Test for randomness
bars_nn <- qgis_run_algorithm(alg="native:nearestneighbouranalysis",
INPUT=results.proj,
OUTPUT_HTML_FILE=file.path(tempdir(), "barsnn.html"),
load_output = TRUE)
## Ignoring unknown input 'load_output'
## Running cmd.exe /c call \
## "C://Program Files//QGIS 3.16.16//bin//qgis_process-qgis-ltr.bat" --json \
## run "native:nearestneighbouranalysis" \
## "--INPUT=C:\Users\anami\AppData\Local\Temp\Rtmp0y7DLJ\file54b8132f451b\file54b8b46515d.gpkg" \
## "--OUTPUT_HTML_FILE=C:\Users\anami\AppData\Local\Temp\Rtmp0y7DLJ/barsnn.html"
## <string>:1: DeprecationWarning: setapi() is deprecated
bars_nn
## <Result of `qgis_run_algorithm("native:nearestneighbouranalysis", ...)`>
## List of 6
## $ EXPECTED_MD : num 1885
## $ NN_INDEX : num 0.963
## $ OBSERVED_MD : num 1815
## $ OUTPUT_HTML_FILE: 'qgis_outputHtml' chr "C:\\Users\\anami\\AppData\\Local\\Temp\\Rtmp0y7DLJ/barsnn.html"
## $ POINT_COUNT : num 25
## $ Z_SCORE : num -0.352
Analysis summary:
Since the nearest neighborhood index value is 0.963 ( less than 1), it exhibits a clustering pattern. Bars are clustered over the geographical area. Also, the Z value is -0.352. As the value is less than -1.96, there is no complete spatial randomness. However, the negative value of Z indicates that a clustered pattern exists.