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.
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:
## [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.