Overview

The AZ House Committee is debating a bill that would require all new wind farm projects to be no closer than 12 miles to property zoned for residential use.

Last year, USA TODAY documented how counties across the US are blocking wind and solar projects, slowing the renewable energy transition.

To understand the impact of the proposed Arizona bill, we mapped out how much land would be left and who owns it.

Here’s what we found:

If 12-mile buffers from residential buildings are applied in Arizona, they would cover around 90% of the state’s land. Roughly half of the land left is in the hands of the National Parks Service (including the Grand Canyon) and in Native American reservations.

Note: No state agency keeps a statewide zoning map. This data is kept as a patchwork at the city and county level so we used FEMA’s USA structures database. This catalogs all buildings greater than 450 square feet and can be filtered by “residential” structures. We calculated the buffers from the center of the buildings due to computer processing limits. At the scales we’re working at, it should give us a good proxy of the impacted areas.

Scroll below for full details on how we got this, summary tables, and interactive maps.

Methods

FEMA data can be downloaded as a geodatabase (gdb) file by state here. Let’s import it in R.

# path to data
gdb_path <- "inputs/Deliverable20230502AZ/AZ_Structures.gdb"

# 2.7M features, 28 fields, crs = WGS 84
st_layers(gdb_path)

# import data
data <- st_read(gdb_path, layer = "AZ_Structures")

There are 2.7 million structures in the state, most of which are residential buildings. These include single/multi family homes, manufactured homes, nursing homes, and dorms.

# summarise by occupancy type
st_drop_geometry(data) %>%
  group_by(OCC_CLS) %>%
  summarise(count = n(),
            pct = round(n()/nrow(data), 3)) %>%
  arrange(-pct) %>%
  adorn_totals()
##           OCC_CLS   count   pct
##       Residential 2368182 0.877
##        Commercial  109725 0.041
##      Unclassified   89560 0.033
##        Government   35291 0.013
##        Industrial   30297 0.011
##         Education   19389 0.007
##       Agriculture   15392 0.006
##          Assembly   17066 0.006
##  Utility and Misc   16889 0.006
##             Total 2701791 1.000
# summarise residential buildings by primary occupancy types 
st_drop_geometry(data) %>%
  filter(OCC_CLS == "Residential") %>%
  group_by(PRIM_OCC) %>%
  summarise(count = n(),
            pct = round(n()/2368182, 3)) %>%
  arrange(-pct) %>%
  adorn_totals()
##                 PRIM_OCC   count   pct
##   Single Family Dwelling 1889127 0.798
##        Manufactured Home  305892 0.129
##  Multi - Family Dwelling  143803 0.061
##             Unclassified   16795 0.007
##             Nursing Home    5206 0.002
##        Temporary Lodging    4597 0.002
##  Institutional Dormitory    2762 0.001
##                    Total 2368182 1.000

Let’s make a copy of the data filtered by the residential buildings. Then we need to re-project the geometries from degrees to meters for accurate calculations.

Getting the buffer around the building perimeters takes too much computing power, so we’ll calculate centroids and then buffer around those.

# first filter data for residential structures
data2 <- filter(data, OCC_CLS == "Residential") %>% select()

# re-project from degrees to meters
data3 <- st_transform(data2, crs = st_crs("ESRI:102010"))

# get centroids and then buffer around those
centroids <- st_centroid(data3$Shape)
buffers <- st_buffer(centroids, dist = (12 * 1609.34)) # meters to miles

buffer_merged <- st_union(buffers) %>% st_sf()

With that area ready, we can get the state’s boundary using the “tigris” package. We’ll use that to clip the buffer and then get the difference.

# now get AZ's boundary shapefile
az <- states(cb = TRUE) %>% 
  filter(STUSPS == "AZ") %>%
  select() %>%
  st_transform(st_crs(buffer_merged)) # set the crs the same as buffer

# clip the buffer edges
buffer_int <- st_intersection(az, buffer_merged)

# now get the difference
buffer_diff <- st_difference(az, buffer_int)

# map full buffer (crs 4326 makes the map level)
ggplot() +
  geom_sf(data = st_transform(az, crs = 4326), color = "black", fill = "orange", alpha = 1, size = 1.2) +
  geom_sf(data = st_transform(buffer_merged, crs = 4326), color = "blue", fill = "blue", alpha = 0.5, size = 1.2) +
  theme_void() +
  labs(title = "12-mile buffer from residential buildings, AZ")

Now we can calculate exactly how much land is left after the buffer.

# square meters to square miles conversion factor
m2_mi2 <- 2589988.11

# state total area (compare to Census Bureau: https://data.census.gov/profile/Arizona?g=040XX00US04)
area_az <- as.numeric(st_area(az))/m2_mi2

# buffer area
area_buffer_int <- as.numeric(st_area(buffer_int))/m2_mi2

# remaining area after buffer removal
area_buffer_diff <- as.numeric(st_area(buffer_diff))/m2_mi2

# pct area of state in buffer
area_buffer_int_pct <- round(area_buffer_int/area_az, 4)

The 12-mile buffer from residential buildings would cover 89.17% of Arizona. Now for some additional context, we can see what type of land is left, which might reduce the available area even further.

The Arizona State Land Department has a land ownership map which we can extract using the ArcGIS server.

# define the ArcGIS REST API endpoint
url <- "http://gis.azland.gov/arcgis/rest/services/Administrative_Boundaries/PublicLandOwnership/MapServer/0/query"

# query parameters
query <- list(where = "1=1", # select all features
              outFields = "*", # retrieve all fields
              f = "geojson",
              returnGeometry = "true", # ensure geometry is included
              geometryType = "esriGeometryPolygon",
              spatialRel = "esriSpatialRelIntersects",
              returnIdsOnly = "false",
              returnCountOnly = "false",
              returnDistinctValues = "false",
              returnTrueCurves = "false")

# fetch data
response <- GET(url, query = query)

# convert response to text and parse JSON into sf
geojson_text <- content(response, as = "text", encoding = "UTF-8")

land_types <- st_read(geojson_text, quiet = TRUE) %>%
  st_make_valid() %>%
  st_transform(st_crs(buffer_diff)) # consistent crs for calculations

rm(query, response, geojson_text, url)

The data includes the following land categories:

levels(as.factor(land_types$category))
##  [1] "BLM"                  "BR"                   "County"              
##  [4] "IBWC"                 "Indian Lands"         "Local or State Parks"
##  [7] "Military"             "NPS"                  "Other"               
## [10] "Private"              "State"                "State Wildlife Area" 
## [13] "USFS"                 "USFWS"

Let’s calculate the which of those intersect the buffered area and what percentage they make from the state’s total land.

# first combine by land category
land_types_merged <- land_types %>%
  group_by(category) %>%
  summarize(geometry = st_union(geometry))

# then get the intersection between buffer_diff and land_types_merged
land_merged_int <- st_intersection(buffer_diff, land_types_merged) %>%
  # for some reason the intersection creates geometry collections
  # extract those polygons and merge them into multipolygons
  # this saves some hassle down the line
  st_collection_extract("POLYGON") %>%
  st_cast("MULTIPOLYGON") %>%
  group_by(category) %>%
  summarize(geometry = st_union(geometry))

# calculate how much categories make up from the state's total land
land_merged_int$int_mi2 <- round(as.numeric(st_area(land_merged_int$geometry))/m2_mi2)

land_merged_int$int_pct <- round(land_merged_int$int_mi2/area_az, 4)

# print a table
land_merged_int <- st_drop_geometry(land_merged_int) %>%
  arrange(-int_mi2)

datatable(land_merged_int, extensions = 'Buttons', options = list(
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf')))

Tribal land and National Park Service make up roughly half of what’s left after the buffer.

We can drill down one more level to get extra details from land ownership.

# then get the intersection between buffer_diff and land_types
land_int <- st_intersection(buffer_diff, land_types) %>%
  st_collection_extract("POLYGON") %>%
  st_cast("MULTIPOLYGON") %>%
  group_by(category, desc_) %>%
  summarize(geometry = st_union(geometry))

# calculate how much categories make up from the state's total land
land_int$int_mi2 <- round(as.numeric(st_area(land_int$geometry))/m2_mi2)

land_int$int_pct <- round(land_int$int_mi2/area_az, 4)

# print a table
land_int <- land_int %>%
  select(category, desc_, int_mi2, int_pct) %>%
  arrange(-int_mi2)

datatable(st_drop_geometry(land_int), extensions = 'Buttons', options = list(
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf')))

Let’s map the area left after the buffer one more time but overlaying all ownership categories that intersect. Click the map for more details.