options(Ncores = 12)
library(mapview)
library(sf)
library(censusxy)
library(dplyr)
library(tmap)

Create point layer from data using latitude and longitude and plot points

addr <- read.csv("C:/Users/chrys/Documents/GitHub/DEM7093/data/drinking_places_dt.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, layer.name = "Drinking Establishments")
#mapshot(mv1, file = paste0(getwd(), "/map1.png"))
mv1

Mean feature - average of bar coordinates

Mean Center of the points is: c(2121292.51040046, 13701736.8982003)

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 = "purple")+
  tm_shape(mean_feature)+
  tm_dots(col = "green", size = .2)

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 = "purple")+
  tm_shape(chull)+
  tm_polygons(col = "grey", alpha = .5)

Alternative Convex Hull Plot

mapview(chull, layer.name = "Convex Hull")+
  mapview(results, col.regions = "green", layer.name = "Drinking Establishments")

Spatial join

library(tidycensus)

#load census tract data
sa_acs<-get_acs(geography = "tract",
                state = "TX",
                county = "Bexar", 
                year = 2019,
                variables = "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 outline
spjoin<-st_join(results.proj, sa_acs2)
#head(spjoin)

tmap_mode("view")
## tmap mode set to interactive viewing
tm_basemap("OpenStreetMap.Mapnik")+
  tm_shape(spjoin, is.master = T)+
  tm_dots("totpop", size = .1, title = "Population/Bar")+
  tm_shape(sa_acs2)+
  tm_polygons(alpha = .1)

Count points in polygons

sa_acs2$nbar<- lengths(st_intersects(sa_acs2, results.proj))

sa_acs2$bar_pc <- 1000*(sa_acs2$nbar/sa_acs2$totpop)

tmap_mode("plot")
## tmap mode set to plotting
m<-tm_basemap("OpenStreetMap.Mapnik")+
  tm_shape(sa_acs2)+
  tm_polygons("bar_pc", title = "Bars / 1,000 People")+
  tm_shape(spjoin, is.master = T)+
  tm_dots( size = .07)

m


Alternative Density Map

mapview(sa_acs2, zcol="bar_pc", layer.name = "Bars/1000 Pop")+
  mapview(results.proj, col.regions = "green", layer.name = "Bars")

Nearest Neighbor analysis

Nearest neighbor analysis to see the distributions of space between features, a measure of clustering or dispersion.

library(spatstat)
bar.pp<-as.ppp(as(results.proj, "Spatial"))

plot(nearest.neighbour(bar.pp))

Average Nearest Neighbor (ANN) 1815/1885 = 0.963
ANN less than 1 indicates clustering
A z-score of -0.352 also indicates clustering and is not significant.

library(qgisprocess)
qgis_configure()

bars_nn<-qgis_run_algorithm(alg="native:nearestneighbouranalysis",
                            INPUT=results.proj,
                            OUTPUT_HTML_FILE=file.path(tempdir(), "barsnn.html"),
                            load_output = TRUE)
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\\chrys\\AppData\\Local\\Temp\\Rtmp0S7yb0/barsnn.html"
##  $ POINT_COUNT     : num 25
##  $ Z_SCORE         : num -0.352