Biodiversity Intactness Index (BII) change in Phoenix, AZ

A geospatial analysis in Python
EDS
Geospatial
Author

Bailey Jørgensen

Published

September 29, 2024

About:

Purpose: In 2021, Maricopa County —home to the Phoenix metropolitan area— was identified as the U.S. county with the most significant increase in developed land since 2001. This rapid urban sprawl has profound implications for biodiversity and the health of surrounding natural ecosystems.

In this notebook, I will investigate the impacts of urban expansion by analyzing a dataset that captures values for the Biodiversity Intactness Index (BII). Apecifically, I will examine changes in BII in the Phoenix county subdivision area between 2017 and 2020, shedding light on how urban growth affects biodiversity over time.

Highlights: 1. Accessing the Microsoft Planetary Computer to Access their Impact Observatory data.

  1. Creating a plot of the Biodiversity Intactness Index within the Phoenix subdivision polygon.

  2. Calculate the percentage of area of the Phoenix area with a BII index of at least .07 in 2017 and 2020.

  3. Plotting the 2020 data, revealing areas that show a loss of biodiversity.

About the data: 1. The first data set is is the Biodiversity Intactness Index (BII) Time Series. Access the io-biodiversity collection from the Microsoft Planetary Computer STAC catalog. I will be using the 2017 and 2020 rasters covering the Phoenix subdivision.

  1. The second data set is the Phoenix Subdivision Shapefile Download the Phoenix subdivision polygon from the Census County Subdivision shapefiles for Arizona. All legal boundaries and names are as of January 1, 2024. The 2024 TIGER/Line Shapefiles were released on September 25, 2024. https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2022&layergroup=County+Subdivisions

Visualization: Map of Phoenix County AZ

In order to better visualize the Phoenix area, I will make a map, and add in a basemap using the Contextily package. This map will give us a good visual of the greater Phoenix area, and will better help us understand our biodiveristy analysis later.

# Set up the figure
fig, ax = plt.subplots()

# Create the axis with plot
phoenix_merc = phoenix.to_crs(epsg=3857)
phoenix_merc.plot(ax=ax, figsize=(11, 10), alpha=0.45, edgecolor="k", label="Phoenix Boundary")

# Add NatGeo basemap from contextily
ctx.add_basemap(ax=ax, source=ctx.providers.Esri.NatGeoWorldMap)

# Add legend
legend_elements = [Line2D([0], [0], color='k', alpha=0.45, lw=3, label='Phoenix County Boundary')]
ax.legend(handles=legend_elements, loc="upper left", fontsize=12)

# Update axes
ax.set_title("Phoenix County, Arizona", fontdict={"fontsize": 20})
ax.set_axis_off()

plt.show()

Biodiversity Intactness Index

Now that we have a good basemap, showing roads and ciites for context, we can work on also plotting our Plantary Computer data, to see biodiversity intactness in this highly urban area.

Remembering the data tables created above, I know that I am interested the “data” asset, which should contain the information that I need to continue my comparative analysis. Since my analysis will involve directly comparing the years 2017 and 2020, I will create variables defining the assets of those years specifically.

# Retrieve 2017 biodiversity data 
phx_2017_asset = phx_2017.assets["data"]
phx_2017_data = rioxr.open_rasterio(phx_2017_asset.href)

# Retrieve 2020 biodiversity data 
phx_2020_asset = phx_2020.assets["data"]
phx_2020_data = rioxr.open_rasterio(phx_2020_asset.href)

phx_2017_data
phx_2020_data
<xarray.DataArray (band: 1, y: 7992, x: 7992)> Size: 255MB
[63872064 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
  * x            (x) float64 64kB -115.4 -115.4 -115.4 ... -108.2 -108.2 -108.2
  * y            (y) float64 64kB 34.74 34.74 34.74 34.74 ... 27.57 27.57 27.57
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

Since there is only 1 band, we can remove that dimension for a more easily managed dataset.

phx_2017_data = phx_2017_data.squeeze().drop_vars('band')
phx_2020_data = phx_2020_data.squeeze().drop_vars('band')

Now, I need to clip the Planetary Computer data to match the Phoenix polygon I created earlier. I will make sure that my Coordinate Reference Systems match, before clipping and plotting.

# Match CRSs
phoenix = phoenix.to_crs(phx_2017_data.rio.crs)

# Assert check to ensure it worked
assert phoenix.crs == phx_2017_data.rio.crs

# Clip 2017 BII raster to the Phoenix polygon
phx_clip_2017 = phx_2017_data.rio.clip(phoenix["geometry"])

# And again for 2020
phx_clip_2020 =  phx_2020_data.rio.clip(phoenix["geometry"])

# Let's see how that looks for 2017
phx_clip_2017.plot()

Calculating the percentage of area of the Phoenix subdivision with a BII of at least 0.75 in 2017 and 2020.

In order to do this calculation, I need to do a bit of “map algebra”.

  1. I need to find the total number of pixels in the Phoenix subdivision polygon, for both years.

  2. Then, I need to find the numbers of pixels with the BII at 0.75 or higher, again for both years.

  3. Finally, I can use those pixel numbers I get to calculate the desired percentage.

# Calculate total area for 2017
total_2017 = phx_clip_2017.count().item() 
# Calculate total area for 2020
total_2020 = phx_clip_2020.count().item()

# Calculate BII % for 2017, Make the data binary for easy math
phx_bii_2017 = (phx_clip_2017 >= 0.75).astype(int)
pixels_2017 = phx_bii_2017.sum().item()

# Calculate BII % for 2020
phx_bii_2020 = (phx_clip_2020 >= 0.75).astype(int)
pixels_2020 = phx_bii_2020.sum().item()



# Calcualte the percentage area for 2017 data:
pct_bii_2017 =  (pixels_2017 / total_2017) * 100

# Calculate the percentage area for 2020 data:
pct_bii_2020 = (pixels_2020 / total_2020) * 100

# Print output
pct_text = "The percentage of area in Phoenix County with a BII over 0.75 in"
print(pct_text, "2017 is: ", round(pct_bii_2017, 2), "%")

print(pct_text, "2020 is: ", round(pct_bii_2020, 2), "%")
The percentage of area in Phoenix County with a BII over 0.75 in 2017 is:  7.13 %
The percentage of area in Phoenix County with a BII over 0.75 in 2020 is:  6.49 %

Visualization: Biodiversity Loss

Now that we have calculated the differences in biodiversity density between the two years, we can visualize this difference with a map.

# Find areas where there was a loss in BII from 2017 to 2020
phx_bii_diff = phx_bii_2017 - phx_bii_2020

# Make a mask - like in 223!
phx_mask = phx_bii_diff == 1

# Convert to int type
phx_mask = phx_mask.astype(int)

# Make a color gradient
color_id = ["none", "red"]

# Make a color map object
cmap = plt.cm.colors.ListedColormap(color_id)

# Plot!
phx_mask.plot(cmap = cmap, add_colorbar = False)

Now that we have applied a color map visualization to our biodiversity loss, we can overlay this with our Planetary Computer raster, and our Phoenix subdivision polygon, to create a final visualization that shows us the area, in red, where Phoenix county experienced significant biodiversity losee.

Final Visualization:

# Initialize plot
fig, ax = plt.subplots(figsize = (8, 7))

# Remove axes
ax.axis("off")

# Begin with clipped raster layer
phx_clip_2020.plot(ax = ax,
                  cbar_kwargs = {"location": "bottom",
                                "label": "BII in 2020"})
# Add mask layer
phx_mask.plot(ax = ax,
             cmap = cmap,
             add_colorbar = False)
phx_mask_patch = mpatches.Patch(color = "red",
                      label = "Area with 75% Biodiversity Loss between 2017 and 2020")
# Add polygon layer
phoenix.plot(ax = ax, 
            color = "none",
            edgecolor = "black",
            linewidth = 2)

# update legend, add in BII loss patch
ax.legend(handles = [phx_mask_patch],
          frameon = False,
          bbox_to_anchor = (0.85, -0.07))
ax.set_title("Biodiversity Intactness Index (BII) and Biodiversity Loss in Phoenix, AZ")

plt.show()

Conclusion

Python and Jupyter Notebooks, as we can see, are powerful tools that allow us to conduct high level analyis with large amounts of spatial data.

In this analysis, we can use these tools to visualize the sharp decrease in biodiversity in this area, perhaps due to urbanization and human impact. Such an analysis can be used to further understand the overall biodiversity crisis, and can help inform decision makers about the needs and priorities of conservation efforts.

Citations and Data Access:

C. Galaz García, EDS 220 - Working with Environmental Datasets, Course Notes. 2024. [Online]. Available: https://meds-eds-220.github.io/MEDS-eds-220-course/book/preface.html

Microsoft Planetary Computer, STAC Catalog. Biodiversity Intactness (‘io-biodiversity’). [Dataset]. https://planetarycomputer.microsoft.com/dataset/io-biodiversity Accessed 2 December 2024.

United States Census Bureau. (2022). Arizona County Subdivision 2022 TIGER/Line Shapefiles. [Data File]. U.S. Census Bureau, Geography Division. https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2022&layergroup=County+Subdivisions Accessed 2 December 2024.

Citation

BibTeX citation:
@online{jørgensen2024,
  author = {Jørgensen, Bailey},
  title = {Biodiversity {Intactness} {Index} {(BII)} Change in
    {Phoenix,} {AZ}},
  date = {2024-09-29},
  url = {https://jorb1.github.io/posts/2024-09-19-BII/},
  langid = {en}
}
For attribution, please cite this work as:
Jørgensen, Bailey. 2024. “Biodiversity Intactness Index (BII) Change in Phoenix, AZ.” September 29, 2024. https://jorb1.github.io/posts/2024-09-19-BII/.