GIS Assignment 6
Jyoti Nepal, MSW
07 March, 2022
knitr::opts_chunk$set(echo = TRUE)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(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)library(readr)
addr<-read_csv("C:/Users/spara/OneDrive/Desktop/New folder/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...##
## 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.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"))
mv1library(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)
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 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
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 viewingtm_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()algs<-qgis_algorithms()
algs[grepl(pattern = "density", x = algs$algorithm ),]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
## Argument type: distance
## Acceptable values:
## - A numeric value
## RADIUS_FIELD: Radius from field
## Argument type: field
## Acceptable values:
## - The name of an existing field
## - ; delimited list of existing field names
## PIXEL_SIZE: Output raster size
## Argument type: number
## Acceptable values:
## - A numeric value
## WEIGHT_FIELD: Weight from field
## Argument type: field
## Acceptable values:
## - The name of an existing field
## - ; delimited list of existing field names
## KERNEL: Kernel shape
## 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)
## Argument type: number
## Acceptable values:
## - A numeric value
## OUTPUT_VALUE: Output value scaling
## 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
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 cmd.exe /c call \
## "C:/Program Files/QGIS 3.22.4/bin/qgis_process-qgis-ltr.bat" --json run \
## "qgis:heatmapkerneldensityestimation" \
## "--INPUT=C:\Users\spara\AppData\Local\Temp\RtmpaCpuA8\file1ec443c73865\file1ec45ffa4737.gpkg" \
## "--RADIUS=2000" "--PIXEL_SIZE=100" "--KERNEL=0" "--OUTPUT_VALUE=0" \
## "--OUTPUT=C:/Users/spara/OneDrive/Desktop/wicdenst.TIF"## <string>:1: DeprecationWarning: setapi() is deprecatedlibrary(raster)## Loading required package: sp##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectlibrary(RColorBrewer)
result<- qgis_as_raster(wic_dens)
projection(result)<-crs(results.proj)
tmap_mode("view")## tmap mode set to interactive viewingtm_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 outlinespjoin<-st_join(results.proj, sa_acs2)
#head(spjoin)
tmap_mode("view")## tmap mode set to interactive viewingtm_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 plottingm<-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 polygonswic_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
map1Nearest 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-scoreAs the neighborhood index is less than 1, the pattern is clustered.