xarray backend: MUR SST¶
The MUR SST dataset has daily records for sea surface temperature and ice cover fraction. There is a netcdf file for each record.
To run the titiler-cmr service locally you can fire up the docker network with this command:
docker compose up
Requirements¶
To run some of the chunks in this notebook you will need to install a few packages:
earthaccessfoliumhttpxxarray
!pip install folium httpx earthaccess xarray
In [15]:
Copied!
import json
from datetime import datetime, timezone
import earthaccess
import httpx
import xarray as xr
from folium import GeoJson, Map, TileLayer
# titiler_endpoint = "http://localhost:8081" # docker network endpoint
titiler_endpoint = "https://dev-titiler-cmr.delta-backend.com" # deployed endpoint
import json
from datetime import datetime, timezone
import earthaccess
import httpx
import xarray as xr
from folium import GeoJson, Map, TileLayer
# titiler_endpoint = "http://localhost:8081" # docker network endpoint
titiler_endpoint = "https://dev-titiler-cmr.delta-backend.com" # deployed endpoint
Identify the dataset¶
You can find the MUR SST dataset using the earthaccess.search_datasets function.
In [16]:
Copied!
datasets = earthaccess.search_datasets(doi="10.5067/GHGMR-4FJ04")
ds = datasets[0]
concept_id = ds["meta"]["concept-id"]
print("Concept-Id: ", concept_id)
print("Abstract: ", ds["umm"]["Abstract"])
datasets = earthaccess.search_datasets(doi="10.5067/GHGMR-4FJ04")
ds = datasets[0]
concept_id = ds["meta"]["concept-id"]
print("Concept-Id: ", concept_id)
print("Abstract: ", ds["umm"]["Abstract"])
Concept-Id: C1996881146-POCLOUD Abstract: A Group for High Resolution Sea Surface Temperature (GHRSST) Level 4 sea surface temperature analysis produced as a retrospective dataset (four day latency) and near-real-time dataset (one day latency) at the JPL Physical Oceanography DAAC using wavelets as basis functions in an optimal interpolation approach on a global 0.01 degree grid. The version 4 Multiscale Ultrahigh Resolution (MUR) L4 analysis is based upon nighttime GHRSST L2P skin and subskin SST observations from several instruments including the NASA Advanced Microwave Scanning Radiometer-EOS (AMSR-E), the JAXA Advanced Microwave Scanning Radiometer 2 on GCOM-W1, the Moderate Resolution Imaging Spectroradiometers (MODIS) on the NASA Aqua and Terra platforms, the US Navy microwave WindSat radiometer, the Advanced Very High Resolution Radiometer (AVHRR) on several NOAA satellites, and in situ SST observations from the NOAA iQuam project. The ice concentration data are from the archives at the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI SAF) High Latitude Processing Center and are also used for an improved SST parameterization for the high-latitudes. The dataset also contains additional variables for some granules including a SST anomaly derived from a MUR climatology and the temporal distance to the nearest IR measurement for each pixel.This dataset is funded by the NASA MEaSUREs program ( http://earthdata.nasa.gov/our-community/community-data-system-programs/measures-projects ), and created by a team led by Dr. Toshio M. Chin from JPL. It adheres to the GHRSST Data Processing Specification (GDS) version 2 format specifications. Use the file global metadata "history:" attribute to determine if a granule is near-realtime or retrospective.
Examine a granule¶
Each granule contains a single day record for the entire globe and has a single data file.
In [17]:
Copied!
results = earthaccess.search_data(
count=1,
concept_id=concept_id,
temporal=("2024-10-12", "2024-10-13"),
)
print("Granules:")
print(results)
print()
print("Example of NetCDF URL: ")
for link in results[0].data_links(access="external"):
print(link)
results = earthaccess.search_data(
count=1,
concept_id=concept_id,
temporal=("2024-10-12", "2024-10-13"),
)
print("Granules:")
print(results)
print()
print("Example of NetCDF URL: ")
for link in results[0].data_links(access="external"):
print(link)
Granules:
[Collection: {'Version': '4.1', 'ShortName': 'MUR-JPL-L4-GLOB-v4.1'}
Spatial coverage: {'HorizontalSpatialDomain': {'Geometry': {'BoundingRectangles': [{'WestBoundingCoordinate': -180, 'SouthBoundingCoordinate': -90, 'EastBoundingCoordinate': 180, 'NorthBoundingCoordinate': 90}]}}}
Temporal coverage: {'RangeDateTime': {'EndingDateTime': '2024-10-12T21:00:00.000Z', 'BeginningDateTime': '2024-10-11T21:00:00.000Z'}}
Size(MB): 707.340648651123
Data: ['https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20241012090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc']]
Example of NetCDF URL:
https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-protected/MUR-JPL-L4-GLOB-v4.1/20241012090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1.nc
Explore the available variables¶
The NetCDF file can be opened with xarray using the h5netcdf engine. When running outside of AWS region us-west-2 you will need to access the data using "external" https links (rather than "direct" s3 links). Those links will require authentication which is handled by earthaccess as long as you have your Earthdata credentials stored in the ~/.netrc file!
In [21]:
Copied!
fs = earthaccess.get_fsspec_https_session()
ds = xr.open_dataset(fs.open(results[0].data_links(access="external")[0]), engine="h5netcdf")
print("Data Variables:")
for var in ds.data_vars:
print(str(var))
display(ds)
fs = earthaccess.get_fsspec_https_session()
ds = xr.open_dataset(fs.open(results[0].data_links(access="external")[0]), engine="h5netcdf")
print("Data Variables:")
for var in ds.data_vars:
print(str(var))
display(ds)
Data Variables: analysed_sst analysis_error mask sea_ice_fraction dt_1km_data sst_anomaly
<xarray.Dataset> Size: 29GB
Dimensions: (time: 1, lat: 17999, lon: 36000)
Coordinates:
* time (time) datetime64[ns] 8B 2024-10-12T09:00:00
* lat (lat) float32 72kB -89.99 -89.98 -89.97 ... 89.98 89.99
* lon (lon) float32 144kB -180.0 -180.0 -180.0 ... 180.0 180.0
Data variables:
analysed_sst (time, lat, lon) float64 5GB ...
analysis_error (time, lat, lon) float64 5GB ...
mask (time, lat, lon) float32 3GB ...
sea_ice_fraction (time, lat, lon) float64 5GB ...
dt_1km_data (time, lat, lon) timedelta64[ns] 5GB ...
sst_anomaly (time, lat, lon) float64 5GB ...
Attributes: (12/47)
Conventions: CF-1.7
title: Daily MUR SST, Final product
summary: A merged, multi-sensor L4 Foundation SST anal...
references: http://podaac.jpl.nasa.gov/Multi-scale_Ultra-...
institution: Jet Propulsion Laboratory
history: created at nominal 4-day latency; replaced nr...
... ...
project: NASA Making Earth Science Data Records for Us...
publisher_name: GHRSST Project Office
publisher_url: http://www.ghrsst.org
publisher_email: [email protected]
processing_level: L4
cdm_data_type: gridDefine a query for titiler-cmr¶
To use titiler-cmr's endpoints for a NetCDF dataset like this we need to define a date range for the CMR query and a variable to analyze.
In [11]:
Copied!
variable = "sea_ice_fraction"
datetime_ = datetime(2024, 10, 10, tzinfo=timezone.utc).isoformat()
variable = "sea_ice_fraction"
datetime_ = datetime(2024, 10, 10, tzinfo=timezone.utc).isoformat()
Display tiles in an interactive map¶
The /tilejson.json endpoint will provide a parameterized xyz tile URL that can be added to an interactive map.
In [12]:
Copied!
r = httpx.get(
f"{titiler_endpoint}/WebMercatorQuad/tilejson.json",
params = (
("concept_id", concept_id),
# Datetime in form of `start_date/end_date`
("datetime", datetime_),
# titiler-cmr can work with both Zarr and COG dataset
# but we need to tell the endpoints in advance which backend
# to use
("backend", "xarray"),
("variable", variable),
# We need to set min/max zoom because we don't want to use lowerzoom level (e.g 0)
# which will results in useless large scale query
("minzoom", 2),
("maxzoom", 13),
("rescale", "0,1"),
("colormap_name", "blues_r"),
)
).json()
print(r)
r = httpx.get(
f"{titiler_endpoint}/WebMercatorQuad/tilejson.json",
params = (
("concept_id", concept_id),
# Datetime in form of `start_date/end_date`
("datetime", datetime_),
# titiler-cmr can work with both Zarr and COG dataset
# but we need to tell the endpoints in advance which backend
# to use
("backend", "xarray"),
("variable", variable),
# We need to set min/max zoom because we don't want to use lowerzoom level (e.g 0)
# which will results in useless large scale query
("minzoom", 2),
("maxzoom", 13),
("rescale", "0,1"),
("colormap_name", "blues_r"),
)
).json()
print(r)
{'tilejson': '2.2.0', 'version': '1.0.0', 'scheme': 'xyz', 'tiles': ['https://dev-titiler-cmr.delta-backend.com/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?concept_id=C1996881146-POCLOUD&datetime=2024-10-10T00%3A00%3A00%2B00%3A00&backend=xarray&variable=sea_ice_fraction&rescale=0%2C1&colormap_name=blues_r'], 'minzoom': 2, 'maxzoom': 13, 'bounds': [-180.0, -90.0, 180.0, 90.0], 'center': [0.0, 0.0, 2]}
In [13]:
Copied!
bounds = r["bounds"]
m = Map(
location=(70, -40),
zoom_start=3
)
TileLayer(
tiles=r["tiles"][0],
opacity=1,
attr="NASA",
).add_to(m)
m
bounds = r["bounds"]
m = Map(
location=(70, -40),
zoom_start=3
)
TileLayer(
tiles=r["tiles"][0],
opacity=1,
attr="NASA",
).add_to(m)
m
Out[13]:
Make this Notebook Trusted to load map: File -> Trust Notebook
GeoJSON Statistics¶
The /statistics endpoint can be used to get summary statistics for a geojson Feature or FeatureCollection.
In [14]:
Copied!
geojson_dict = {
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"coordinates": [
[
[
-20.79973248834736,
83.55979308678764
],
[
-20.79973248834736,
75.0115425216471
],
[
14.483337068956956,
75.0115425216471
],
[
14.483337068956956,
83.55979308678764
],
[
-20.79973248834736,
83.55979308678764
]
]
],
"type": "Polygon"
}
}
]
}
r = httpx.post(
f"{titiler_endpoint}/statistics",
params=(
("concept_id", concept_id),
# Datetime in form of `start_date/end_date`
("datetime", datetime_),
# titiler-cmr can work with both Zarr and COG dataset
# but we need to tell the endpoints in advance which backend
# to use
("backend", "xarray"),
("variable", variable),
),
json=geojson_dict,
timeout=60,
).json()
print(json.dumps(r, indent=2))
geojson_dict = {
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"coordinates": [
[
[
-20.79973248834736,
83.55979308678764
],
[
-20.79973248834736,
75.0115425216471
],
[
14.483337068956956,
75.0115425216471
],
[
14.483337068956956,
83.55979308678764
],
[
-20.79973248834736,
83.55979308678764
]
]
],
"type": "Polygon"
}
}
]
}
r = httpx.post(
f"{titiler_endpoint}/statistics",
params=(
("concept_id", concept_id),
# Datetime in form of `start_date/end_date`
("datetime", datetime_),
# titiler-cmr can work with both Zarr and COG dataset
# but we need to tell the endpoints in advance which backend
# to use
("backend", "xarray"),
("variable", variable),
),
json=geojson_dict,
timeout=60,
).json()
print(json.dumps(r, indent=2))
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
-20.79973248834736,
83.55979308678764
],
[
-20.79973248834736,
75.0115425216471
],
[
14.483337068956956,
75.0115425216471
],
[
14.483337068956956,
83.55979308678764
],
[
-20.79973248834736,
83.55979308678764
]
]
]
},
"properties": {
"statistics": {
"sea_ice_fraction": {
"min": 0.3,
"max": 0.99,
"mean": 0.845157064600111,
"count": 1725290.875,
"sum": 1458141.771496357,
"std": 0.1559272507275522,
"median": 0.9,
"majority": 0.9500000000000001,
"minority": 0.36,
"unique": 70.0,
"histogram": [
[
34892,
39574,
38696,
37867,
44348,
72817,
110580,
200188,
472678,
675707
],
[
0.3,
0.369,
0.43799999999999994,
0.5069999999999999,
0.576,
0.645,
0.714,
0.7829999999999999,
0.8519999999999999,
0.9209999999999998,
0.99
]
],
"valid_percent": 57.18,
"masked_pixels": 1293477.0,
"valid_pixels": 1727347.0,
"percentile_2": 0.36,
"percentile_98": 0.99
}
}
}
}
]
}
In [ ]:
Copied!