Comparing NLDAS-2 and NLDAS-3 Precipitation Forcing Data

Author

Rishi Anand, Siddharth Chaudhary

Published

August 29, 2024

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'

About the Data

NLDAS is a widely used land modeling environment that generates estimates of land surface fluxes and states such as soil moisture, snow, and streamflow. These estimates are critical for drought and flood monitoring, water availability and water resource management, climate assessments, and other uses. NLDAS-3 is the next generation version of NLDAS-2, and offers significant improvements such as improved spatial resolution (12.5km to 1km), expanded domain (CONUS to North and Central America), reduced data latency (3.5 days to near real-time), and assimilation of NASA remote sensing data, among others. (see Earthdata VEDA Data Story). Please note that the NLDAS-3 precipitation data provided here is a sample dataset still in development, and will not be the final NLDAS-3 product when it is released.

This notebook is intended to visualize and compare the NLDAS-2 and sample NLDAS-3 monthly-averaged precipitation forcing.

Approach

  1. Query metadata from the VEDA STAC API
  2. Get tiles representing NLDAS-2 and NLDAS-3 from the VEDA Raster API
  3. Display the data side-by-side for comparison
import requests
import folium
import folium.plugins

from pystac.client import Client

Querying the STAC API

The metadata and derived data products can be accessed without credentials via STAC API and Raster API respectively.

STAC_API_URL = "https://openveda.cloud/api/stac"
RASTER_API_URL = "https://openveda.cloud/api/raster"

Next we will define our collections of interest. You can discover available collections via the STAC Browser: http://openveda.cloud

collection_id3 = "nldas3"
collection_id2 = "nldas2"

Now we can use those collection ids along with a particular datetime to search the STAC for items that match our criteria.

search = Client.open(STAC_API_URL).search(
    collections=[collection_id3, collection_id2],
    datetime="2021-02-01",
)

item_collection = search.item_collection()
print(f"Found {len(item_collection)} items")
Found 2 items

Lets take a closer look at each of these items:

item_collection[0]
item_collection[1]

Notice that each of these items has an asset called cog_default. That is the one we are interested in. That asset includes statistics about the raster band:

nldas3_item = item_collection[0]
nldas3_band_statistics = nldas3_item.assets["cog_default"].ext.raster.bands[0].statistics
nldas3_band_statistics.properties
{'mean': 30.96841488575868,
 'stddev': 41.008411449510426,
 'maximum': 428.0843200683594,
 'minimum': -0.022795898839831352,
 'valid_percent': 47.08281935307018}

Get tiles with the Raster API

We’ll use those statistics to make a Raster API query for the NLDAS-3 tile:

response = requests.get(
    f"{RASTER_API_URL}/collections/{collection_id3}/items/{nldas3_item.id}/WebMercatorQuad/tilejson.json?"
    "&assets=cog_default"
    "&color_formula=gamma+r+1.05&colormap_name=rdbu_r"
    f"&rescale={nldas3_band_statistics.minimum},{nldas3_band_statistics.maximum}",
)

# If anything goes wrong with this request, output error contents
assert response.ok, response.text
    
tiles3 = response.json()
tiles3
{'tilejson': '2.2.0',
 'version': '1.0.0',
 'scheme': 'xyz',
 'tiles': ['https://openveda.cloud/api/raster/collections/nldas3/items/nldas3_LIS_HIST_202102/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?assets=cog_default&color_formula=gamma+r+1.05&colormap_name=rdbu_r&rescale=-0.022795898839831352%2C428.0843200683594'],
 'minzoom': 0,
 'maxzoom': 24,
 'bounds': [-168.98000671647551,
  7.019999961559591,
  -51.93999145246978,
  72.06000091582078],
 'center': [-110.45999908447266, 39.540000438690186, 0]}

Let’s do the same for the nldas2 item

nldas2_item = item_collection[1]
nldas2_band_statistics = nldas2_item.assets["cog_default"].ext.raster.bands[0].statistics

response = requests.get(
    f"{RASTER_API_URL}/collections/{collection_id2}/items/{nldas2_item.id}/WebMercatorQuad/tilejson.json?"
    "&assets=cog_default"
    "&color_formula=gamma+r+1.05&colormap_name=rdbu_r"
    f"&rescale={nldas2_band_statistics.minimum},{nldas2_band_statistics.maximum}",
)

# If anything goes wrong with this request, output error contents
assert response.ok, response.text

tiles2 = response.json()
tiles2
{'tilejson': '2.2.0',
 'version': '1.0.0',
 'scheme': 'xyz',
 'tiles': ['https://openveda.cloud/api/raster/collections/nldas2/items/nldas2_LIS_HIST_202102/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?assets=cog_default&color_formula=gamma+r+1.05&colormap_name=rdbu_r&rescale=-4.244315147399902%2C574.3709106445312'],
 'minzoom': 0,
 'maxzoom': 24,
 'bounds': [-168.98000671647551,
  7.019999961559591,
  -51.93999145246978,
  72.06000091582078],
 'center': [-110.45999908447266, 39.540000438690186, 0]}

Visualize NLDAS-2 and NLDAS-3

Create and display a DualMap to visualize NLDAS-2 and NLDAS-3 data side-by-side. Notice the differences in extent and resolution between the two versions.

map_layer2 = folium.TileLayer(
    tiles=tiles2["tiles"][0],
    attr="VEDA",
)
map_layer3 = folium.TileLayer(
    tiles=tiles3["tiles"][0],
    attr="VEDA",
)
m = folium.plugins.DualMap(location=[39,-110], zoom_start=3,)
map_layer2.add_to(m.m1)
map_layer3.add_to(m.m2)
m
Make this Notebook Trusted to load map: File -> Trust Notebook