Homework 6
Mahmuda Sultana
3/3/2022
library(mapview)
library(sf)## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUElibrary(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, unionlibrary(readr)
Drinking_places <- read_csv("C:/Users/shahi/OneDrive - University of Texas at San Antonio/Desktop/Spring'22/Dem 5093 (GIS)/Drinking_places.csv")## New names:
## * Source -> Source...1
## * Source -> Source...227
## * `` -> ...228## Rows: 25 Columns: 228## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (76): Source...1, Date, Obsolescence Date, Business Name, Legal Name, P...
## dbl (40): Physical Address Number, Physical ZIP, Physical ZIP 4, Location E...
## lgl (112): Physical Post Direction, Corporate Employee Size, Mailing Post Di...##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.Drinking_places <- Drinking_places%>%
dplyr::select(names(Drinking_places)[c(4,6,14,41,42)])
results <- st_as_sf(Drinking_places, 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"))
mv1 webshot::install_phantomjs(force = TRUE)## phantomjs has been installed to C:\Users\shahi\AppData\Roaming\PhantomJSlibrary(rJava)
library(tmap)
library(tmaptools)
library(OpenStreetMap)bg<- read_osm(results.proj, ext=1.1)
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots()tmap_mode("plot")## tmap mode set to plottinglibrary(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.library(ggspatial)
library(leaflet)
leaflet() %>%
addTiles()%>%
addMarkers(lng = Drinking_places$Longitude, lat= Drinking_places$Latitude)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)
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj, size = .2)+
tm_dots(col = "green")+
tm_shape(mean_feature)+
tm_dots(col = "red", size = .2)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)
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col = "green", size = .2)+
tm_shape(mean_feature)+
tm_dots(col = "red", size=.2)+
tm_shape(median_feature)+
tm_dots(col = "blue", size = .2)Buffer points
Drinking_places_buff<- st_buffer(results.proj, dist = 1000)
#tmap_mode("plot")
tm_basemap("OpenStreetMap.Mapnik" )+
tm_shape(results.proj, is.master = T)+
tm_dots(col = "green")+
tm_shape(Drinking_places_buff)+
tm_polygons(col = "red", alpha = .1)Convex hull plot
chull <- st_convex_hull(st_union(results))
tmap_mode("view")## tmap mode set to interactive viewingtm_basemap("OpenStreetMap.Mapnik")+
tm_shape(results.proj)+
tm_dots(col = "green")+
tm_shape(chull)+
tm_polygons(col = "grey", alpha = .5)kernel density
R can do kernel density maps, but using simple features it’s kind of complicated. Here is one way to do this using the SpatialKDE package:
library(SpatialKDE)
grid_drinking <- results.proj %>%
create_grid_rectangular(cell_size = 1000, side_offset = 2000)
kde <- results.proj%>%
kde(band_width = 2000, kernel= "quartic", grid = grid_drinking)
tm_shape(kde)+
tm_polygons(col = "kde_value", palette= "viridis", title = "Density Estimate")+
tm_shape(results.proj)+
tm_dots()Spatial join
A spatial join can combine attributes of one layer with another layer. Here I combine census variables with the Drinking places.
library(tidycensus)
#load census tract data
sa_acs<-get_acs(geography = "tract",
state="TX",
county = "Bexar",
year = 2019,
variables=c( "DP05_0001E") ,
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#rename variables and filter missing cases
sa_acs2<-sa_acs%>%
mutate(totpop= DP05_0001E)%>%
dplyr::select(GEOID, totpop)
sa_acs2<-st_transform(sa_acs2, crs = 2278)
sa_trol<-st_cast(sa_acs2, "MULTILINESTRING") #make outlinespjoin<-st_join(results.proj, sa_acs2)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$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 plottingm<-tm_basemap("OpenStreetMap.Mapnik")+
tm_shape(sa_acs2)+
tm_polygons("bars_pc")+
tm_shape(spjoin, is.master = T)+
tm_dots( size = .2)
m
Nearest Neighbor analysis
Nearest neighbor analysis is used to test if a pattern of points is distrubuted randomly or clustered.
library(spatstat)
drinking.pp<-as.ppp(as(results.proj, "Spatial"))
plot(nearest.neighbour(drinking.pp))
options(qgisprocess.path = "C://Program Files//QGIS 3.16.15//bin//qgis_process-qgis-ltr.bat" )
library(qgisprocess)## Using 'qgis_process' at 'C://Program Files//QGIS 3.16.15//bin//qgis_process-qgis-ltr.bat'.
## QGIS version: 3.16.15-Hannover
## Configuration loaded from 'C:\Users\shahi\AppData\Local/R-qgisprocess/R-qgisprocess/Cache/cache-0.0.0.9000.rds'
## Run `qgis_configure()` for details.
## Using JSON for output serializationqgis_configure()## Trying getOption('qgisprocess.path'): 'C://Program Files//QGIS 3.16.15//bin//qgis_process-qgis-ltr.bat'## Success!## QGIS version: 3.16.15-Hannover## Saving configuration to 'C:\Users\shahi\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 serializationHere I use QGIS to get a single test of the 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.15//bin//qgis_process-qgis-ltr.bat" --json \
## run "native:nearestneighbouranalysis" \
## "--INPUT=C:\Users\shahi\AppData\Local\Temp\RtmpsdUL2O\file11fc42c447b\file11fc18216a4e.gpkg" \
## "--OUTPUT_HTML_FILE=C:\Users\shahi\AppData\Local\Temp\RtmpsdUL2O/barsnn.html"## <string>:1: DeprecationWarning: setapi() is deprecatedbars_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\\shahi\\AppData\\Local\\Temp\\RtmpsdUL2O/barsnn.html"
## $ POINT_COUNT : num 25
## $ Z_SCORE : num -0.352##Report the numeric summary of the analysis:
Here, the nearest neighbour index is 0.963, which is less than 1. It indicates clustering pattern. Z-sore is -0.352, this negative z-score also indicates clustering but it’s insignificant as it’s less than 1.96.