Land Use/Land Cover

This page was generated from notebooks/nlcd.ipynb. Interactive online version: Binder badge

Land Use/Land Cover

import networkx as nx
import osmnx as ox
import pygeohydro as gh
from pynhd import NLDI

Land cover, imperviousness, and canopy data can be retrieved from the NLCD database. First, we use PyNHD to get the contributing watershed geometry of three NWIS stations with the ID of USGS-01031450, USGS-01031500, and USGS-01031510:

geometry = NLDI().get_basins(["01031450", "01031500", "01031510"])

We can now use the nlcd_bygeom and nlcd_bycoords functions to get the NLCD data.

Let’s start by nlcd_bygeom. This function has two positional arguments for passing the target geometries or points of interests and target resolution in meters. Note that, if a single geometry is passed and the geometry is not in EPSG:4326 CRS, geo_crs argument should be given as well. The second argument is the target resolution of the data in meters. The NLCD database is multi-resolution and based on the target resolution, the source data are resampled on the server side.

You should be mindful of the resolution since higher resolutions require more memory so if your requests requires more memory than the available memory on your system the code is likely to crash. You can either increase the resolution or divide your region of interest into smaller regions.

Moreover, the MRLC GeoServer has a limit of about 8 million pixels per request but PyGeoHydro takes care of the domain decomposition under-the-hood and divides the request to smaller requests then merges them. So the only bottleneck for requests is the amount of available memory on your system.

Both nlcd_bygeom and nlcd_bycoords functions can request for four layers from the MRLC web service; imperviousness, land use/land cover, impervious descriptors, and tree canopy. Since NLCD is released every couple of years, you can specify the target year via the years argument. Layers that are not included in this argument are ignored. By default, years is {'impervious': [2019], 'cover': [2019], 'canopy': [2019], 'descriptor': [2019]}.

Furthermore, we can specify the region of interest as well via the region argument. Valid values are L48 (for CONUS), HI (for Hawaii), AK (for Alaska), and PR (for Puerto Rico). By default, region is set to L48.

Let’s get the cover and impervious descriptor data at a 100 m resolution for all three stations

desc = gh.nlcd_bygeom(geometry, 100, years={"descriptor": 2019})

This function returns a dict where its keys are the indices of the input GeoDataFrame.

cmap, norm, levels = gh.plot.descriptor_legends()
ax = desc["01031500"].descriptor_2019.plot(
    size=5, cmap=cmap, levels=levels, cbar_kwargs={"ticks": levels[:-1]}
ax.axes.set_title("Urban Descriptor 2019")
ax.figure.savefig("_static/descriptor.png", bbox_inches="tight", facecolor="w")

Now let’s get the land cover data:

lulc = gh.nlcd_bygeom(geometry, 100, years={"cover": [2016, 2019]})

Additionally, PyGeoHydro provides a function for getting the official legends of the cover data. Let’s plot the data using this legends.

cmap, norm, levels = gh.plot.cover_legends()
ax = lulc["01031500"].cover_2019.plot(
    size=5, cmap=cmap, levels=levels, cbar_kwargs={"ticks": levels[:-1]}
ax.axes.set_title("Land Use/Land Cover 2019")
ax.figure.savefig("_static/lulc.png", bbox_inches="tight", facecolor="w")

Moreover, we can get the statistics of the cover data for each class or category as well:

stats = gh.cover_statistics(lulc["01031500"].cover_2019)
{'classes': {'Open Water': 2.232661383991947,
  'Developed, Open Space': 2.471414190950623,
  'Developed, Low Intensity': 0.7872389851069871,
  'Developed, Medium Intensity': 0.32522003974911595,
  'Developed, High Intensity': 0.06323722995121699,
  'Barren Land (Rock/Sand/Clay)': 0.03355444854554371,
  'Deciduous Forest': 20.463051389928502,
  'Evergreen Forest': 21.975582685904552,
  'Mixed Forest': 36.216864982061274,
  'Shrub-Forest': 3.4251348630720386,
  'Herbaceous-Forest': 1.4570374003045712,
  'Shrub/Scrub': 0.17551557700745943,
  'Grassland/Herbaceous': 0.16131946416126786,
  'Pasture/Hay': 1.546085744521591,
  'Cultivated Crops': 0.06581834137779728,
  'Woody Wetlands': 7.993702088119144,
  'Emergent Herbaceous Wetlands': 0.6065611852463672},
 'categories': {'Background': 0.0,
  'Unclassified': 0.0,
  'Water': 2.232661383991947,
  'Developed': 3.6471104457579435,
  'Barren': 0.03355444854554371,
  'Forest': 83.53767132127093,
  'Shrubland': 0.17551557700745943,
  'Herbaceous': 0.16131946416126786,
  'Planted/Cultivated': 1.6119040858993885,
  'Wetlands': 8.600263273365512}}

Now, let’s see nlcd_bycoords in action. The coordinates must be a list of (longitude, latitude) coordinates. Let’s use osmnx package to get a street network:

G = ox.graph_from_place("Piedmont, California, USA", network_type="drive")

Now, we can get land cover and tree canopy for each node based on their coordinates and then plot the results.

x, y = nx.get_node_attributes(G, "x").values(), nx.get_node_attributes(G, "y").values()
lulc = gh.nlcd_bycoords(list(zip(x, y)), years={"cover": [2019], "canopy": [2016]})
nx.set_node_attributes(G, dict(zip(G.nodes(), lulc.cover_2019)), "cover_2019")
nx.set_node_attributes(G, dict(zip(G.nodes(), lulc.canopy_2016)), "canopy_2016")
geometry canopy_2016 cover_2019
0 POINT (-122.24760 37.82625) 11.0 22
1 POINT (-122.24716 37.82420) 0.0 23
2 POINT (-122.24607 37.82491) 0.0 23
3 POINT (-122.24532 37.82541) 0.0 22
4 POINT (-122.24447 37.82595) 5.0 23
nc = ox.plot.get_node_colors_by_attr(G, "canopy_2016", cmap="viridis_r")
fig, ax = ox.plot_graph(