library(mapview)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(censusxy)
## Warning: package 'censusxy' was built under R version 4.0.4
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(sf)
library(readxl)
## Warning: package 'readxl' was built under R version 4.0.4
WIC <- read_excel("C:/Users/adolp/Desktop/GIS/WS_Stores.xls")
## New names:
## * Source -> Source...1
## * Source -> Source...227
results <- st_as_sf(WIC, coords=c("Longitude", "Latitude"), crs=4269,agr="constant")
results.proj<-st_transform(results,
crs = 2278)
mapview(results.proj)
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)
mapview(mean_feature, col.regions="red")+mapview( results)
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)
mapview(median_feature, col.regions="green")+
mapview(mean_feature, col.regions="red")+
mapview( results)
Storebuff<- st_buffer(results.proj, dist = 2500)
mapview(Storebuff)+mapview(results.proj, col.regions="green")
chull <- st_convex_hull(st_union(results))
## although coordinates are longitude/latitude, st_union assumes that they are planar
mapview(chull)+
mapview(results, col.regions = "green")
remotes::install_github("paleolimbot/qgisprocess")
## Skipping install of 'qgisprocess' from a github remote, the SHA1 (f065b10c) has not changed since last install.
## Use `force = TRUE` to force installation
library(qgisprocess)
## Using 'qgis_process' at 'C:/Program Files/QGIS 3.18/bin/qgis_process-qgis.bat'.
## Run `qgis_configure()` for details.
qgis_configure()
## getOption('qgisprocess.path') was not found.
## Sys.getenv('R_QGISPROCESS_PATH') was not found.
## Trying 'qgis_process' on PATH
## Error in rethrow_call(c_processx_exec, command, c(command, args), pty, : Command 'qgis_process' not found @win/processx.c:982 (processx_exec)
## Found 1 QGIS installation containing 'qgis_process':
## C:/Program Files/QGIS 3.18/bin/qgis_process-qgis.bat
## Trying command 'C:/Program Files/QGIS 3.18/bin/qgis_process-qgis.bat'
## Success!
algs<-qgis_algorithms()
algs[grepl(pattern = "density", x = algs$algorithm ),]
## # A tibble: 8 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 grass7 GRASS grass7:r.li.ed~ r.li.edgedensity r.li.edgedensity
## 2 grass7 GRASS grass7:r.li.ed~ r.li.edgedensit~ r.li.edgedensity.as~
## 3 grass7 GRASS grass7:r.li.pa~ r.li.patchdensi~ r.li.patchdensity
## 4 grass7 GRASS grass7:r.li.pa~ r.li.patchdensi~ r.li.patchdensity.a~
## 5 native QGIS (native c~ native:lineden~ linedensity Line density
## 6 qgis QGIS qgis:heatmapke~ heatmapkernelde~ Heatmap (Kernel Den~
## 7 saga SAGA saga:fragmenta~ fragmentationcl~ Fragmentation class~
## 8 saga SAGA saga:kernelden~ kerneldensityes~ Kernel density esti~
qgis_show_help("qgis:heatmapkerneldensityestimation")
## Heatmap (Kernel Density Estimation) (qgis:heatmapkerneldensityestimation)
##
## ----------------
## Description
## ----------------
##
## ----------------
## 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>
## Heatmap
Store_dens<-qgis_run_algorithm(algorithm ="qgis:heatmapkerneldensityestimation",
INPUT=results.proj,
RADIUS = 5000,
PIXEL_SIZE = 100,
KERNEL = 0,
OUTPUT=file.path(getwd(), "Storedenst.TIF"),
load_output = TRUE)
## Running "C:/Program Files/QGIS 3.18/bin/qgis_process-qgis.bat" run \
## "qgis:heatmapkerneldensityestimation" \
## "--INPUT=C:\Users\adolp\AppData\Local\Temp\RtmpSGHIkg\file338026d17642\file3380475c2ae1.gpkg" \
## "--RADIUS=5000" "--PIXEL_SIZE=100" "--KERNEL=0" "--OUTPUT_VALUE=0" \
## "--OUTPUT=C:/Users/adolp/Desktop/GIS/Storedenst.TIF"
##
## ----------------
## Inputs
## ----------------
##
## INPUT: C:\Users\adolp\AppData\Local\Temp\RtmpSGHIkg\file338026d17642\file3380475c2ae1.gpkg
## KERNEL: 0
## OUTPUT: C:/Users/adolp/Desktop/GIS/Storedenst.TIF
## OUTPUT_VALUE: 0
## PIXEL_SIZE: 100
## RADIUS: 5000
##
##
## 0...10...20...30...40...50...60...70...80...90...
## ----------------
## Results
## ----------------
##
## OUTPUT: C:/Users/adolp/Desktop/GIS/Storedenst.TIF
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
library(RColorBrewer)
result<- qgis_as_raster(Store_dens)
projection(result)<-crs(results.proj)
mapview(result)+mapview(results.proj)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in CRS definition
library(tidycensus)
library(dplyr)
sa_acs<-get_acs(geography = "tract",
state="TX",
county = "Bexar",
year = 2019,
variables=c( "DP05_0001E", "DP03_0009P", "DP03_0062E", "DP03_0119PE",
"DP05_0001E","DP02_0009PE","DP02_0008PE","DP02_0040E","DP02_0038E",
"DP02_0066PE","DP02_0067PE","DP02_0080PE","DP02_0092PE",
"DP03_0005PE","DP03_0028PE","DP03_0062E","DP03_0099PE","DP03_0101PE",
"DP03_0119PE","DP04_0046PE","DP05_0072PE","DP05_0073PE",
"DP05_0066PE", "DP05_0072PE", "DP02_0113PE") ,
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
## Using the ACS Data Profile
sa_acs2<-sa_acs%>%
mutate(totpop= DP05_0001E, pwhite=DP05_0072PE,
pblack=DP05_0073PE , phisp=DP05_0066PE,
phsormore=DP02_0066PE,punemp=DP03_0009PE, medhhinc=DP03_0062E,
ppov=DP03_0119PE)%>%
dplyr::select(GEOID, totpop, pblack, pwhite, phisp, punemp, medhhinc, ppov)
sa_acs2<-st_transform(sa_acs2, crs = 2278)
sa_trol<-st_cast(sa_acs2, "MULTILINESTRING")
spjoin<-st_join(results.proj, sa_acs2)
## New names:
## * Source...227 -> Source...225
head(spjoin)
## Simple feature collection with 6 features and 233 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 2111381 ymin: 13697280 xmax: 2128787 ymax: 13726580
## projected CRS: NAD83 / Texas South Central (ftUS)
## # A tibble: 6 x 234
## Source...1 Date `Obsolescence Da~ `Business Name` `Legal Name`
## <chr> <chr> <chr> <chr> <chr>
## 1 AtoZDatabas~ 03/24/2~ 09/24/2021 H-E-B LP HEB Grocery Com~
## 2 AtoZDatabas~ 03/24/2~ 09/24/2021 Surlean Foods <NA>
## 3 AtoZDatabas~ 03/24/2~ 09/24/2021 Dean Foods <NA>
## 4 AtoZDatabas~ 03/24/2~ 09/24/2021 Super Target <NA>
## 5 AtoZDatabas~ 03/24/2~ 09/24/2021 H-E-B LP <NA>
## 6 AtoZDatabas~ 03/24/2~ 09/24/2021 Marios & Brother Cor~ Mario's Grocery
## # ... with 229 more variables: Physical Address <chr>,
## # Physical Address Number <chr>, Physical Pre Direction <chr>,
## # Physical Address Name <chr>, Physical Address Suffix <chr>,
## # Physical Post Direction <lgl>, Physical City <chr>, Physical State <chr>,
## # Physical ZIP <chr>, Physical ZIP 4 <chr>, Key Executive Name <chr>,
## # First Name <chr>, Middle Initial <chr>, Last Name <chr>, Title <chr>,
## # Gender <chr>, Location Employee Size <chr>, Corporate Employee Size <chr>,
## # Revenue / Yr <chr>, Mailing Address <chr>, Mailing Address Number <chr>,
## # Mailing Pre Direction <chr>, Mailing Address Name <chr>,
## # Mailing Address Suffix <chr>, Mailing Post Direction <lgl>,
## # Mailing City <chr>, Mailing State <chr>, Mailing ZIP <chr>,
## # Mailing ZIP 4 <chr>, Phone <chr>, Fax <chr>, Toll-Free <chr>,
## # County Name <chr>, County Population <chr>, Metro Area <chr>, EIN <chr>,
## # Main Line of Business <chr>, Location Type <chr>,
## # Importer or Exporter <lgl>, Manufacturer <chr>, Primary SIC <chr>,
## # Primary SIC Description <chr>, SIC02 <chr>, SIC02.Description <chr>,
## # SIC03 <chr>, SIC03.Description <chr>, SIC04 <chr>, SIC04.Description <chr>,
## # SIC05 <chr>, SIC05.Description <chr>, SIC06 <chr>, SIC06.Description <chr>,
## # SIC07 <chr>, SIC07.Description <chr>, SIC08 <chr>, SIC08.Description <chr>,
## # SIC09 <chr>, SIC09.Description <chr>, SIC10 <chr>, SIC10.Description <chr>,
## # NAICS 1 <chr>, NAICS 1 Description <chr>, NAICS 2 <chr>,
## # NAICS 2 Description <chr>, NAICS 3 <chr>, NAICS 3 Description <chr>,
## # NAICS 4 <chr>, NAICS 4 Description <chr>, NAICS 5 <chr>,
## # NAICS 5 Description <chr>, NAICS 6 <chr>, NAICS 6 Description <chr>,
## # NAICS 7 <chr>, NAICS 7 Description <chr>, NAICS 8 <chr>,
## # NAICS 8 Description <chr>, NAICS 9 <chr>, NAICS 9 Description <chr>,
## # NAICS 10 <lgl>, NAICS 10 Description <lgl>, Non-Profit <chr>,
## # Number of PCs <chr>, Public / Private <chr>, Small Business <chr>,
## # Square Footage <chr>, Website <chr>, Women Owned <chr>,
## # Year Established <chr>, Ticker Symbol <lgl>, Stock Exchange <lgl>,
## # Fortune 1000 Ranking <lgl>, Credit Score <chr>, 2020 Revenue/Yr <chr>,
## # 2019 Revenue/Yr <chr>, 2018 Revenue/Yr <chr>, 2017 Revenue/Yr <chr>,
## # 2019 % Sales Growth <chr>, 2018 % Sales Growth <chr>,
## # 2017 % Sales Growth <chr>, 2020 Employees <chr>, ...
mapview(spjoin["punemp"])+mapview(sa_trol)
sa_acs2$nwic<- lengths(st_intersects(sa_acs2, results.proj))
mapview(sa_acs2, zcol="nwic")+
mapview(results.proj, col.regions = "green")
algs[grepl(pattern = "voronoi", x = algs$algorithm ),]
## # A tibble: 3 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 grass7 GRASS grass7:v.voronoi v.voronoi v.voronoi
## 2 grass7 GRASS grass7:v.voronoi.sk~ v.voronoi.skele~ v.voronoi.skele~
## 3 qgis QGIS qgis:voronoipolygons voronoipolygons Voronoi polygons
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
Store_von<-qgis_run_algorithm(alg="qgis:voronoipolygons",
INPUT=results.proj,
OUTPUT=file.path(tempdir(), "Storevon.shp"),
load_output = TRUE)
Store_von<-sf::read_sf(qgis_output(Store_von, "OUTPUT"))
mapview(Store_von, alpha=.75)+
mapview(results.proj, col.regions="green")
library(spatstat)
Store.pp<-as.ppp(as(results.proj, "Spatial"))
plot(nearest.neighbour(Store.pp))
algs[grepl(pattern = "nearest", x = algs$algorithm ),]
## # A tibble: 10 x 5
## provider provider_title algorithm algorithm_id algorithm_title
## <chr> <chr> <chr> <chr> <chr>
## 1 gdal GDAL gdal:gridinver~ gridinversedist~ Grid (IDW with nea~
## 2 gdal GDAL gdal:gridneare~ gridnearestneig~ Grid (Nearest neig~
## 3 native QGIS (native c~ native:angleto~ angletonearest Align points to fe~
## 4 native QGIS (native c~ native:joinbyn~ joinbynearest Join attributes by~
## 5 native QGIS (native c~ native:nearest~ nearestneighbou~ Nearest neighbour ~
## 6 qgis QGIS qgis:distancet~ distancetoneare~ Distance to neares~
## 7 qgis QGIS qgis:distancet~ distancetoneare~ Distance to neares~
## 8 qgis QGIS qgis:knearestc~ knearestconcave~ Concave hull (k-ne~
## 9 saga SAGA saga:knearestn~ knearestneighbo~ K-nearest neighbou~
## 10 saga SAGA saga:nearestne~ nearestneighbour Nearest neighbour
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
Store_nn<-qgis_run_algorithm(alg="native:nearestneighbouranalysis",
INPUT=results.proj,
OUTPUT_HTML_FILE=file.path(tempdir(), "Storenn.html"),
load_output = TRUE)
## Ignoring unknown input 'load_output'
## Running "C:/Program Files/QGIS 3.18/bin/qgis_process-qgis.bat" run \
## "native:nearestneighbouranalysis" \
## "--INPUT=C:\Users\adolp\AppData\Local\Temp\RtmpSGHIkg\file338026d17642\file33801aeb36c3.gpkg" \
## "--OUTPUT_HTML_FILE=C:\Users\adolp\AppData\Local\Temp\RtmpSGHIkg/Storenn.html"
##
## ----------------
## Inputs
## ----------------
##
## INPUT: C:\Users\adolp\AppData\Local\Temp\RtmpSGHIkg\file338026d17642\file33801aeb36c3.gpkg
## OUTPUT_HTML_FILE: C:\Users\adolp\AppData\Local\Temp\RtmpSGHIkg/Storenn.html
##
##
## 0...10...20...30...40...50...60...70...80...90...100 - done.
##
## ----------------
## Results
## ----------------
##
## EXPECTED_MD: 1654.757340237235
## NN_INDEX: 0.8046634820351503
## OBSERVED_MD: 1331.5228033185174
## OUTPUT_HTML_FILE: C:\Users\adolp\AppData\Local\Temp\RtmpSGHIkg/Storenn.html
## POINT_COUNT: 59
## Z_SCORE: -2.8703861806878748