This page was generated from examples/nlcd.ipynb. Interactive online version:

NLCD

from pynhd import NLDI

import pygeohydro as gh

import warnings

warnings.filterwarnings("ignore")

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

geometry = NLDI().get_basins("01031500").geometry[0]

We can now use the nlcd function to get the data. This function has one positional argumnet for passing the geometry which could be a Polygon or a boundig box. If CRS of the geometry is not EPSG:4326, it should be passed via geo_crs argument. The second argumnet 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 avialable 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.

The nlcd function can request for three layers from the MRLC web service; imperviousness, land use/land cover, and tree canopy. Since NLCD is released every couple of years, you can specify the target year via the years argument. By default it is 2016. You can set any of the years to None to turn off that layer.

Let’s get the cover data at a 100 m resolution.

lulc = gh.nlcd(geometry, 100, years={"impervious": None, "cover": 2016, "canopy": None})

Additionaly, 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.cover.plot(
    size=5, cmap=cmap, levels=levels, cbar_kwargs={"ticks": levels[:-1]}
)
ax.figure.savefig("../_static/lulc.png", bbox_inches="tight", dpi=100)

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

stats = gh.cover_statistics(lulc.cover)
stats

{'classes': {'Open Water': 2.3100947267893557,
  'Perennial Ice/Snow': 0.0,
  'Developed, Open Space': 2.839222569238314,
  'Developed, Low Intensity': 0.5923650724001756,
  'Developed, Medium Intensity': 0.153576129881527,
  'Developed, High Intensity': 0.04387889425186486,
  'Barren Land (Rock/Sand/Clay)': 0.07098056423095786,
  'Unconsolidated Shore': 0.0,
  'Deciduous Forest': 18.86534341687531,
  'Evergreen Forest': 21.05154479518881,
  'Mixed Forest': 33.87837802957954,
  'Shrub': 7.756239836873757,
  'Herbaceous': 2.181039155460341,
  'Dwarf Scrub': 0.0,
  'Shrub/Scrub': 0.29424670263015257,
  'Grassland/Herbaceous': 0.1780966884340397,
  'Sedge/Herbaceous': 0.0,
  'Lichens': 0.0,
  'Moss': 0.0,
  'Pasture/Hay': 1.4079962831995456,
  'Cultivated Crops': 0.06452778566450713,
  'Woody Wetlands': 7.650414268383966,
  'Emergent Herbaceous Wetlands': 0.6620550809178432},
 'categories': {'Unclassified': 0.0,
  'Water': 2.3100947267893557,
  'Developed': 3.6290426657718813,
  'Barren': 0.07098056423095786,
  'Forest': 83.73254523397775,
  'Shrubland': 0.29424670263015257,
  'Herbaceous': 0.1780966884340397,
  'Planted/Cultivated': 1.4725240688640528,
  'Wetlands': 8.31246934930181}}