Impact to Weld County from Increasing Buffer to 2000 ft (in R)

Given the recent news of Colorado’s regulators supporting a 2,000 ft buffer around existing homes, I thought it would be a good time to take a look at the impact using the sf package in R.

Follow My Blog

Packages

Here are the packages we will need:

1. sf – Simple features (geometry) in R. Basically treats them like a dataframe.
2. dplyr – Data manipulation, part of the tidyverse
3. leaflet – Interactive maps in R
4. tigris – Pull in county shapefiles from US Census
5. downloader – Helpful downloading package

---

Data

All I really need is a GIS layer of home addresses, as well as a shapefile of Weld County.

Luckily, the address data is available openly. Of course, there is probably a large portion of these addresses that are not actually living dwellings, but it’s probably still close, so just keep that in mind when reviewing the data.

#If Not Installed
install.packages('downloader')

#Load downloader
library(downloader)
url1 <- 'https://opendata.arcgis.com/datasets/1483cb91e097462383929f6aecb61649_0.zip'
downloader::download(url1, destfile = 'dataset.zip')
unzip('dataset.zip')
unlink('dataset.zip')

library(dplyr)
library(tigris)
library(sf)

#Read our shapefile in
shp1 <- read_sf('Address_Points_open_data.shp')

#Grab Weld county shapefile
weld <- tigris::counties(state = 'CO')
weld <- st_as_sf(weld)
weld <- weld %>% filter(NAME == 'Weld')

---

2000 Ft Setback

Now we want to transform each sf dataframe to the crs for Colorado North in XY (2231), since we are trying to calculate a specific radius.

Next, I will do some transformations. When you see several nested functions, you work from the inside out.

1. st_buffer -> Creating a 2000 ft buffer zone around each address point.
2. st_union -> Combine all of the “circles” created by buffer into shapefiles.
3. st_intersection -> Find everywhere that the new shapefile and Weld County intersect. Because I select Weld first, it condenses it down to a single row multi-polygon object, which is easier and quicker to map in leaflet. Though still not fast.
4. st_transform -> Convert back to WGS 84 (4326)

This one does take a little bit of time to run, so be aware.

weld <- weld %>% st_transform(crs = 2231)

shp1 <- shp1 %>% st_transform(crs = 2231)

tst2 <- st_transform(
  st_intersection(weld, st_union(
    st_buffer(
      shp1, dist = 2000))), 
        4326)

sum(st_area(tst2$geometry))/4046.86

The last part calculates how many net acres are in the impacted zone, or 925,153 acres. We can also plot the impacted zone.

library(leaflet)
leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addProviderTiles(providers$OpenStreetMap) %>%
  addPolygons(data = tst2$geometry,  color = 'black', fillColor = 'green', fillOpacity = 0.4, weight = 0.4)

Visually, everywhere that is impacted is in green below.

Figure 1: Area Impacted by 2000 ft setback – Weld County

---

Original Impacted Area

I still want to compare to the original setback (200 ft), so I’ll just rerun the calculation but change the dist to be 200 ft.

tst2 <- st_transform(
  st_intersection(weld, st_union(
    st_buffer(
      shp1, dist = 200))),
        4326)

sum(st_area(tst2$geometry))/4046.86

The original impacted area was 83,451 acres, so we are seeing ~841,702 new net acres that fall within our new catchment areas. I’d map it but the file is much much larger and takes a while to render, so feel free to try on your own.

---

Summary

Short and simple, but operators are right to be concerned. Most of the good acreage is actually in the Southwestern portion of Weld, and that is where the majority of acreage is impacted.