Visualizing NCEO Aboveground Woody Biomass 2017 prioritization areas

Explores NCEO Aboveground Woody Biomass priority areas in Guinea using zonal statistics.

Run this notebook

You can launch this notebook in VEDA JupyterHub by clicking the link below.

Launch in VEDA JupyterHub (requires access)

Learn more

Inside the Hub

This notebook was written on the VEDA JupyterHub and as such is designed to be run on a jupyterhub which is associated with an AWS IAM role which has been granted permissions to the VEDA data store via its bucket policy. The instance used provided 16GB of RAM.

See (VEDA Analytics JupyterHub Access)[https://nasa-impact.github.io/veda-docs/veda-jh-access.html] for information about how to gain access.

Outside the Hub

The data is in a protected bucket. Please request access by emailng aimee@developmentseed.org or alexandra@developmentseed.org and providing your affiliation, interest in or expected use of the dataset and an AWS IAM role or user Amazon Resource Name (ARN). The team will help you configure the cognito client.

You should then run:

%run -i 'cognito_login.py'

Approach

1. Query STAC API and explore item contents for a given collection
2. Read and access the data
3. Visualize the collection with hvplot
4. Run zonal statistics on collection using rasterstats
5. Visualize resultant zonal statistics on a choropleth map

About the Data

The NCEO Aboveground Woody Biomass 2017 dataset is a map for Africa at 100 m spatial resolution which was developed using a combination of LiDAR, Synthetic Aperture Radar (SAR) and optical based data. Aboveground woody biomass (AGB) plays an key role in the study of the Earth’s carbon cycle and response to climate change. Estimation based on Earth Observation measurements is an effective method for regional scale studies and the results are expressed as dry matter in Mg ha-1.

Important Note: Users of this dataset should keep in mind that the map is a continental-scale dataset, generated using a combination of different remote sensing data types, with a single method for the whole study area produced in 2017. Users, therefore, should understand that accuracy may vary for different regions and vegetation types.

The Case Study - Guinea

Mapping and understanding the spatial distribution of AGB is key to understanding carbon dioxide emissions from tropical deforestation through the loss of woody carbon stocks. The resulting carbon fluxes from these land-use changes and vegetation degradation can have negative impacts on the global carbon cycle. Change analysis between AGB maps overtime can display losses in high biomass forests, due to suspected deforestation and forest degredation.

The forests of southern Guinea are reported to have some of the highest density AGB of any forest in the world and are one of the most threatened ecoregions in Africa. Importantly, this area was also the epicenter of the 2014 Ebola outbreak, which had a lasting impact on the region. There is more and more evidence that human deforestation activities in this area may have accelerated the spread of the deadly virus as a result of increasing human-bat interactions in the region.

In this example we explore the NCEO AGB dataset for 2017, running zonal statistics at the district (administrative 2) level to understand those areas in Guinea that need greatest prioritization for protection and conservation.

Setting up the Environment

To run zonal statistics we’ll need to import a python package called rasterstats into our environment. You can uncomment the following line for installation. This cell needs only needs to be run once.

!pip install rasterstats --quiet

Querying the STAC API

from pystac_client import Client

# Provide STAC API endpoint
STAC_API_URL = "https://openveda.cloud/api/stac/"

# Declare collection of interest - NCEO Biomass
collection = "nceo_africa_2017"

Now let’s check how many total items are available.

search = Client.open(STAC_API_URL).search(collections=[collection])
items = list(search.items())
print(f"Found {len(items)} items")

Found 1 items

This makes sense as there is only one item available: a map for 2017.

# Explore the "cog_default" asset of one item to see what it contains
items[0].assets["cog_default"].to_dict()

{'href': 's3://nasa-maap-data-store/file-staging/nasa-map/nceo-africa-2017/AGB_map_2017v0m_COG.tif',
 'type': 'image/tiff; application=geotiff; profile=cloud-optimized',
 'title': 'Default COG Layer',
 'description': 'Cloud optimized default layer to display on map',
 'raster:bands': [{'scale': 1,
   'nodata': 'inf',
   'offset': 0,
   'sampling': 'area',
   'data_type': 'uint16',
   'histogram': {'max': 429,
    'min': 0,
    'count': 11,
    'buckets': [405348,
     44948,
     18365,
     6377,
     3675,
     3388,
     3785,
     9453,
     13108,
     1186]},
   'statistics': {'mean': 37.58407913145342,
    'stddev': 81.36678677343947,
    'maximum': 429,
    'minimum': 0,
    'valid_percent': 50.42436439336373}}],
 'roles': ['data', 'layer']}

Explore through the item’s assets. We can see from the data’s statistics values that the min and max values for the observed values range from 0 to 429 Mg ha-1.

Reading and accessing the data

Now that we’ve explored the dataset through the STAC API, let’s read and access the dataset itself.

import odc.stac
import rioxarray


da = odc.stac.load(items[:1], chunks={"latitude": "auto", "longitude": "auto"}, dtype="float32")['cog_default'].squeeze()
da

<xarray.DataArray 'cog_default' (latitude: 81024, longitude: 78077)> Size: 25GB
dask.array<getitem, shape=(81024, 78077), dtype=float32, chunksize=(5792, 5792), chunktype=numpy.ndarray>
Coordinates:
  * latitude     (latitude) float64 648kB 37.73 37.73 37.73 ... -35.05 -35.05
  * longitude    (longitude) float64 625kB -18.27 -18.27 -18.27 ... 51.86 51.86
    spatial_ref  int32 4B 4326
    time         datetime64[ns] 8B 2017-01-01
Attributes:
    nodata:   inf

xarray.DataArray

'cog_default'

• latitude: 81024
• longitude: 78077

• dask.array<chunksize=(5792, 5792), meta=np.ndarray>

  Array Chunk Bytes 23.57 GiB 127.97 MiB Shape (81024, 78077) (5792, 5792) Dask graph 196 chunks in 4 graph layers Data type float32 numpy.ndarray

• Coordinates: (4)
  • latitude
    (latitude)
    float64
    37.73 37.73 37.73 ... -35.05 -35.05

    units :

    degrees_north

    resolution :

    -0.0008983152841195214

    crs :

    EPSG:4326

    array([ 37.730589,  37.729691,  37.728793, ..., -35.051813, -35.052712,
           -35.05361 ], shape=(81024,))

  • longitude
    (longitude)
    float64
    -18.27 -18.27 ... 51.86 51.86

    units :

    degrees_east

    resolution :

    0.0008983152841195214

    crs :

    EPSG:4326

    array([-18.27308 , -18.272182, -18.271284, ...,  51.861987,  51.862885,
            51.863784], shape=(78077,))

  • spatial_ref
    ()
    int32
    4326

    spatial_ref :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    crs_wkt :

    GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    grid_mapping_name :

    latitude_longitude

    GeoTransform :

    -18.273529509559306660548827 0.0008983152841195214287 0 37.731038563588171541596239 0 -0.0008983152841195214287

    array(4326, dtype=int32)

  • time
    ()
    datetime64[ns]
    2017-01-01

    array('2017-01-01T00:00:00.000000000', dtype='datetime64[ns]')

• Attributes: (1)

  nodata :

  inf

In this example, we’ll explore the data contained in the NCEO AGB collection and analyze it for each of the districts in Guinea. To do this we will need to import district (administrative level 2) boundary layers from below. We will use the Humanitarian Data Exchange (HDX) site to retrieve subnational administrative boundaries for Guinea. Specifically, we will use the geoBoundaries-GIN-ADM2_simplified.geojson which can be accessed here and read them in directly using geopandas.

import geopandas as gpd

admin2_gdf = gpd.read_file(
    "https://raw.githubusercontent.com/wmgeolab/geoBoundaries/0f0b6f5fb638e7faf115f876da4e77d8f7fa319f/releaseData/gbOpen/GIN/ADM2/geoBoundaries-GIN-ADM2_simplified.geojson"
)

# check the CRS
print(admin2_gdf.crs)

EPSG:4326

Now we can use the bounds of the admin boundaries to clip the data to a box containing Guinea.

biomass = da.rio.clip_box(*admin2_gdf.total_bounds).compute()

import hvplot.xarray


biomass.hvplot(
    x="longitude",
    y="latitude",
    coastline=True,
    rasterize=True,
    cmap="viridis",
    widget_location="bottom",
    frame_width=600,
)

Zonal Statistics

This map we created above is great, but let’s focus on which districts (administrative level 2 boundaries) should be prioritized for forest conservation.

Zonal statistics is an operation that calculates statistics on the cell values of a raster layer (e.g., the NCEO AGB dataset) within the zones (i.e., polygons) of another dataset. It is an analytical tool that can calculate the mean, median, sum, minimum, maximum, or range in each zone. The zonal extent, often polygons, can be in the form of objects like administrative boundaries, water catchment areas, or field boundaries.

import pandas as pd
from rasterstats import zonal_stats


admin2_biomass = pd.DataFrame(
    zonal_stats(
        admin2_gdf,
        biomass.values,
        affine=biomass.rio.transform(),
        nodata=biomass._FillValue,
        band=1,
    ),
    index=admin2_gdf.index,
)
admin2_biomass

	min	max	mean	count
0	0.0	565.0	51.455244	1276418
1	0.0	546.0	46.844499	527243
2	0.0	513.0	48.882245	1107040
3	0.0	345.0	31.289991	41491
4	0.0	494.0	47.810397	116496
5	0.0	422.0	57.551729	530210
6	0.0	548.0	48.579208	317185
7	0.0	438.0	43.736140	1176850
8	0.0	448.0	44.735186	403305
9	0.0	533.0	78.729734	1315384
10	0.0	422.0	46.984165	426274
11	0.0	452.0	50.073605	161605
12	0.0	483.0	52.073439	1145990
13	0.0	563.0	89.607722	429236
14	0.0	503.0	51.309326	1769535
15	0.0	558.0	65.572246	955756
16	0.0	486.0	48.211255	897398
17	0.0	590.0	90.080966	630855
18	0.0	413.0	47.659512	364101
19	0.0	339.0	37.295832	544563
20	0.0	501.0	57.186423	1614318
21	0.0	470.0	46.781563	216374
22	0.0	469.0	57.429324	278976
23	0.0	592.0	71.860556	461576
24	0.0	623.0	123.975206	814867
25	0.0	451.0	45.048732	859201
26	0.0	564.0	74.174900	1042972
27	0.0	406.0	33.335759	1191395
28	0.0	592.0	86.496324	413448
29	0.0	560.0	57.593215	460749
30	0.0	392.0	29.209417	1813359
31	0.0	554.0	57.995929	771380
32	0.0	389.0	49.646580	615127
33	0.0	588.0	131.347109	326027

Now we’ll join the administrative level 2 boundaries to the zonal statistics results, so that we can map the districts on a choropleth map.

concat_df = admin2_gdf.join(admin2_biomass)
concat_df.head()

	OBJECTID	ISO Code	shapeName	Level	shapeID	shapeGroup	shapeType	geometry	min	max	mean	count
0	1	GN-BE	Beyla	ADM2	GIN-ADM2-49546643B63767081	GIN	ADM2	POLYGON ((-8.24559 8.44255, -8.24158 8.45044, ...	0.0	565.0	51.455244	1276418
1	2	GN-BF	Boffa	ADM2	GIN-ADM2-49546643B69790359	GIN	ADM2	MULTIPOLYGON (((-13.77147 9.84445, -13.76994 9...	0.0	546.0	46.844499	527243
2	3	GN-BK	Boke	ADM2	GIN-ADM2-49546643B67680147	GIN	ADM2	MULTIPOLYGON (((-14.57512 10.76872, -14.57633 ...	0.0	513.0	48.882245	1107040
3	4	GN-C	Conakry	ADM2	GIN-ADM2-49546643B26553537	GIN	ADM2	MULTIPOLYGON (((-13.78686 9.46592, -13.79013 9...	0.0	345.0	31.289991	41491
4	5	GN-CO	Coyah	ADM2	GIN-ADM2-49546643B29309121	GIN	ADM2	POLYGON ((-13.49399 9.53945, -13.4805 9.55304,...	0.0	494.0	47.810397	116496

By sorting the results, we can identify those top districts with the highest mean AGB.

concat_df_sorted = concat_df.sort_values(by="mean", ascending=False)
concat_df_sorted.head()

	OBJECTID	ISO Code	shapeName	Level	shapeID	shapeGroup	shapeType	geometry	min	max	mean	count
33	34	GN-YO	Yomou	ADM2	GIN-ADM2-49546643B32761429	GIN	ADM2	POLYGON ((-9.34981 7.75681, -9.34896 7.7535, -...	0.0	588.0	131.347109	326027
24	25	GN-MC	Macenta	ADM2	GIN-ADM2-49546643B91718973	GIN	ADM2	POLYGON ((-8.95774 8.77472, -9.01024 8.79308, ...	0.0	623.0	123.975206	814867
17	18	GN-KS	Kissidougou	ADM2	GIN-ADM2-49546643B39508892	GIN	ADM2	POLYGON ((-10.45426 9.10945, -10.45334 9.08925...	0.0	590.0	90.080966	630855
13	14	GN-GU	Gueckedou	ADM2	GIN-ADM2-49546643B59147082	GIN	ADM2	POLYGON ((-10.59971 9.05848, -10.59402 9.05494...	0.0	563.0	89.607722	429236
28	29	GN-NZ	Nzerekore	ADM2	GIN-ADM2-49546643B97455025	GIN	ADM2	POLYGON ((-8.93454 8.25441, -8.93687 8.25503, ...	0.0	592.0	86.496324	413448

Visualizing the results with a choropleth map

Now, let’s visualize the results!

import hvplot.pandas

# renaming the shapeName to District for improved legend
concat_df.rename(columns={"shapeName": "District"}, inplace=True)

agb = concat_df.hvplot(
    c="mean",
    width=900,
    height=500,
    geo=True,
    hover_cols=["mean", "District"],
    cmap="viridis",
    hover_fill_color="white",
    line_width=1,
    title="Mean Aboveground Woody Biomass per Guinean District (Mg ha-1)",
    tiles="CartoLight",
)

agb

By hovering over the map, we can identify the names and mean AGB per district.

Summary

In this case study we have successfully performed zonal statistics on the NCEO AGB dataset in Guinea and displayed the results on a choropleth map. The results of this analysis can dispaly those districts which contain the greatest average amount of AGB and should be prioritized for forest protection efforts.