Corey Sparks R blog

I post regularly on various R topics, mostly involving data on people

Intro to Point Pattern Analysis using R and QGIS

Written on March 17, 2021

This example shows how to use R and QGIS from within R to perform a series of common point pattern analysis techniques.

library(mapview)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
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
addr<-read.csv(url("https://raw.githubusercontent.com/coreysparks/data/master/wic_west_side.csv"))
addr<-addr[c(6, 12:14)]
names(addr)<-c("street", "city", "st", "zip")
head(addr)
##                  street        city st   zip
## 1 1039 W Hildebrand Ave San Antonio TX 78201
## 2       1102 Delgado St San Antonio TX 78207
## 3      1115 W Martin St San Antonio TX 78207
## 4      1115 W Martin St San Antonio TX 78207
## 5      1115 W Martin St San Antonio TX 78207
## 6      1115 W Martin St San Antonio TX 78207
results<-cxy_geocode(addr,
                     street = "street",
                     city = "city",
                     state ="st",
                     zip = "zip",
                     class="sf",
                     output = "simple")
## 33 rows removed to create an sf object. These were addresses that the geocoder could not match.
results.proj<-st_transform(results,
                           crs = 2278)

OR just use the lat / long information in the data!

addr<-read.csv(url("https://raw.githubusercontent.com/coreysparks/data/master/wic_west_side.csv"))
results <- st_as_sf(addr, 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"))

library(tmap)
library(tmaptools)
library(OpenStreetMap)
bg<-  read_osm(results.proj, ext=1.1)

tmap_mode("plot")
## tmap mode set to plotting
tm_basemap("OpenStreetMap.Mapnik")+
  tm_shape(results.proj)+
  tm_dots()

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)