Assignment6

library(mapview)
library(sf)

## Linking to GEOS 3.8.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

library(readr)
addr<- read_csv("~/Library/Group Containers/UBF8T346G9.OneDriveSyncClientSuite/OneDrive - University of Texas at San Antonio.noindex/OneDrive - University of Texas at San Antonio/Applied Demography/First Year/Spring 2022/Dem 5093- Corey/drinking.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...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

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)

 webshot::install_phantomjs()

## It seems that the version of `phantomjs` installed is greater than or equal to the requested version.To install the requested version or downgrade to another version, use `force = TRUE`.

mv1 <- mapview(results.proj)
mapshot(mv1, file = paste0(getwd(), "/map1.png"))


mv1

library(tmap)
library(tmaptools)
library(OpenStreetMap)
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.

library(ggspatial)
library(leaflet)
leaflet() %>%
  addTiles()%>%
  addMarkers(lng = addr$Longitude, lat= addr$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 viewing

tm_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 viewing

tm_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

grocery_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(grocery_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 viewing

tm_basemap("OpenStreetMap.Mapnik")+
  tm_shape(results.proj)+
  tm_dots(col = "green")+
  tm_shape(chull)+
  tm_polygons(col = "grey", alpha = .5)

kernel density -

library(SpatialKDE)

grid_groc <- results.proj %>%
  create_grid_rectangular(cell_size = 1000, side_offset = 2000)

kde <- results.proj%>%
  kde(band_width = 2000, kernel= "quartic", grid = grid_groc)

tm_shape(kde)+
  tm_polygons(col = "kde_value", palette= "viridis", title = "Density Estimate")+
  tm_shape(results.proj)+
  tm_dots()

Run the algorithm

wic_dens<-qgis_run_algorithm(algorithm ="qgis:heatmapkerneldensityestimation",
         INPUT=results.proj,
         RADIUS = 2000,
         PIXEL_SIZE = 100,
         KERNEL = 0,
         OUTPUT=file.path(getwd(), "wicdenst.TIF"),
         load_output = TRUE)

## Ignoring unknown input 'load_output'

## Argument `RADIUS_FIELD` is unspecified (using QGIS default value).

## Argument `WEIGHT_FIELD` is unspecified (using QGIS default value).

## Argument `DECAY` is unspecified (using QGIS default value).

## Using `OUTPUT_VALUE = "Raw"`

## Running /Applications/QGIS-LTR.app/Contents/MacOS/bin/qgis_process --json run \
##   'qgis:heatmapkerneldensityestimation' \
##   '--INPUT=/var/folders/8t/qvp10nlx35388lr2r8g61tmw0000gn/T//RtmpacFmka/file9043dd327d3/file9044bdb7dac.gpkg' \
##   '--RADIUS=2000' '--PIXEL_SIZE=100' '--KERNEL=0' '--OUTPUT_VALUE=0' \
##   '--OUTPUT=/Users/josephjaiyeola/Library/CloudStorage/OneDrive-UniversityofTexasatSanAntonio/Applied Demography/First Year/Spring 2022/Dem 5093- Corey/Class_GIS/wicdenst.TIF'

## proj_create_from_database: SQLite error on SELECT name, type, coordinate_system_auth_name, coordinate_system_code, datum_auth_name, datum_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM geodetic_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name

## proj_create_from_database: SQLite error on SELECT name, type, coordinate_system_auth_name, coordinate_system_code, datum_auth_name, datum_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM geodetic_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name

## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/5482

## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/5936

## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/3978
## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/2193

## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/5482

## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/5936

## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/3978
## ERROR 1: PROJ: proj_create_from_database: SQLite error on SELECT name, coordinate_system_auth_name, coordinate_system_code, geodetic_crs_auth_name, geodetic_crs_code, conversion_auth_name, conversion_code, area_of_use_auth_name, area_of_use_code, text_definition, deprecated FROM projected_crs WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## ERROR 1: Cannot parse CRS http://www.opengis.net/def/crs/EPSG/0/2193

## ERROR 1: SQLite error on SELECT name, ellipsoid_auth_name, ellipsoid_code, prime_meridian_auth_name, prime_meridian_code, area_of_use_auth_name, area_of_use_code, publication_date, deprecated FROM geodetic_datum WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## Warning 1: Unable to parse srs_id '2278' well-known text 'PROJCS["NAD83 / Texas South Central (ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",27.8333333333333],PARAMETER["central_meridian",-99],PARAMETER["standard_parallel_1",30.2833333333333],PARAMETER["standard_parallel_2",28.3833333333333],PARAMETER["false_easting",1968500],PARAMETER["false_northing",13123333.333],UNIT["US survey foot",0.304800609601219],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2278"]]'
## ERROR 1: SQLite error on SELECT name, ellipsoid_auth_name, ellipsoid_code, prime_meridian_auth_name, prime_meridian_code, area_of_use_auth_name, area_of_use_code, publication_date, deprecated FROM geodetic_datum WHERE auth_name = ? AND code = ?: no such column: area_of_use_auth_name
## Warning 1: Unable to parse srs_id '2278' well-known text 'PROJCS["NAD83 / Texas South Central (ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",27.8333333333333],PARAMETER["central_meridian",-99],PARAMETER["standard_parallel_1",30.2833333333333],PARAMETER["standard_parallel_2",28.3833333333333],PARAMETER["false_easting",1968500],PARAMETER["false_northing",13123333.333],UNIT["US survey foot",0.304800609601219],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2278"]]'

library(raster)

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

library(RColorBrewer)

result<- qgis_as_raster(wic_dens)

projection(result)<-crs(results.proj)

tmap_mode("view")

## tmap mode set to interactive viewing

tm_basemap("OpenStreetMap.Mapnik")+
  tm_shape(result)+
  tm_raster()+
  tm_shape(results.proj, is.master = T)+
  tm_dots(col="red")

Spatial join

library(tidycensus)

#load census tract data
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

#rename variables and filter missing cases
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") #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("ppov", size = .1)+
  tm_shape(sa_acs2)+
  tm_polygons(alpha = .1)

Count points in polygons

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

sa_acs2$groc_pc <- 1000*(sa_acs2$ngroc/sa_acs2$totpop)

tmap_mode("plot")

## tmap mode set to plotting

m<-tm_basemap("OpenStreetMap.Mapnik")+
  tm_shape(sa_acs2)+
  tm_polygons("groc_pc")+
  tm_shape(spjoin, is.master = T)+
  tm_dots( size = .01)

m

Thiessen/Voronoi Polygons

algs[grepl(pattern = "voronoi", x = algs$algorithm ),]

qgis_show_help("qgis:voronoipolygons")

## Voronoi polygons (qgis:voronoipolygons)
## 
## ----------------
## Description
## ----------------
## This algorithm takes a points layer and generates a polygon layer containing the voronoi polygons corresponding to those input points.
## 
## 
## ----------------
## Arguments
## ----------------
## 
## INPUT: Input layer
##  Argument type:  source
##  Acceptable values:
##      - Path to a vector layer
## BUFFER: Buffer region (% of extent)
##  Argument type:  number
##  Acceptable values:
##      - A numeric value
## OUTPUT: Voronoi polygons
##  Argument type:  sink
##  Acceptable values:
##      - Path for new vector layer
## 
## ----------------
## Outputs
## ----------------
## 
## OUTPUT: <outputVector>
##  Voronoi polygons

wic_von<-qgis_run_algorithm(alg="qgis:voronoipolygons",
         INPUT=results.proj,
         OUTPUT=file.path(tempdir(), "wicvon.shp"),
         load_output = TRUE)

wic_von<-sf::read_sf(qgis_output(wic_von, "OUTPUT"))

tmap_mode("view")
map1<-tm_basemap("OpenStreetMap.Mapnik")+
  tm_shape(wic_von)+
  tm_polygons(col="grey", alpha=.4)+
  tm_shape(results.proj)+
  tm_dots( size = .01)
 

# plot map
map1

Nearest Neighbor analysis

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

plot(nearest.neighbour(wic.pp))

algs[grepl(pattern = "nearest", x = algs$algorithm ),]

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