Coral Reef Watch SST anomaly (griddap -> raster)¶

A griddap dataset returns gridded fields as raster NetCDF. Here we pull three days of NOAA Coral Reef Watch 5 km sea-surface-temperature anomaly over a small equatorial-Pacific box, aggregate each day to a GeoTIFF via the pyramids flow, and read the result back — all without importing xarray.

Setup¶

In [ ]:

import tempfile
from pathlib import Path

from earthlens.aggregate import AggregationConfig
from earthlens.core import EarthLens

import tempfile from pathlib import Path from earthlens.aggregate import AggregationConfig from earthlens.core import EarthLens

Request + download¶

aggregate= is accepted because the dataset is griddap (raster). op='auto' reduces this state field by the daily mean.

In [ ]:

out_dir = Path(tempfile.mkdtemp(prefix='earthlens-erddap-crw-'))
tifs = EarthLens(
    data_source='erddap',
    dataset='NOAA_DHW',
    variables=['CRW_SSTANOMALY'],
    start='2023-06-01',
    end='2023-06-03',
    lat_lim=[0.0, 5.0],
    lon_lim=[150.0, 155.0],
    path=out_dir,
).download(aggregate=AggregationConfig(freq='1D'))

for p in tifs:
    print(f'{p.name}  ({p.stat().st_size // 1024} KB)')

out_dir = Path(tempfile.mkdtemp(prefix='earthlens-erddap-crw-')) tifs = EarthLens( data_source='erddap', dataset='NOAA_DHW', variables=['CRW_SSTANOMALY'], start='2023-06-01', end='2023-06-03', lat_lim=[0.0, 5.0], lon_lim=[150.0, 155.0], path=out_dir, ).download(aggregate=AggregationConfig(freq='1D')) for p in tifs: print(f'{p.name} ({p.stat().st_size // 1024} KB)')

Open a day with pyramids and summarise¶

The fill value (a large negative sentinel) is masked to a physical SST-anomaly range before computing statistics.

In [ ]:

from pyramids.dataset import Dataset

ds = Dataset.read_file(str(tifs[0]))
arr = ds.read_array().astype('float32')
valid = arr[(arr > -50) & (arr < 50)]  # drop the NoData sentinel

print(f'EPSG:        {ds.epsg}')
print(f'array shape: {arr.shape}  ({arr.size:,} pixels)')
print(f'valid px:    {valid.size:,}')
print(f'min anomaly: {valid.min():+.2f} degC')
print(f'mean anomaly:{valid.mean():+.2f} degC')
print(f'max anomaly: {valid.max():+.2f} degC')
print(f'pct > +1 degC (warm stress): {(valid > 1).mean() * 100:.1f} %')

from pyramids.dataset import Dataset ds = Dataset.read_file(str(tifs[0])) arr = ds.read_array().astype('float32') valid = arr[(arr > -50) & (arr < 50)] # drop the NoData sentinel print(f'EPSG: {ds.epsg}') print(f'array shape: {arr.shape} ({arr.size:,} pixels)') print(f'valid px: {valid.size:,}') print(f'min anomaly: {valid.min():+.2f} degC') print(f'mean anomaly:{valid.mean():+.2f} degC') print(f'max anomaly: {valid.max():+.2f} degC') print(f'pct > +1 degC (warm stress): {(valid > 1).mean() * 100:.1f} %')

Daily area-mean anomaly across the window¶

In [ ]:

for p in sorted(tifs):
    a = Dataset.read_file(str(p)).read_array().astype('float32')
    v = a[(a > -50) & (a < 50)]
    day = p.stem.split('_')[-1]
    print(f'{day}: mean anomaly {v.mean():+.3f} degC')

for p in sorted(tifs): a = Dataset.read_file(str(p)).read_array().astype('float32') v = a[(a > -50) & (a < 50)] day = p.stem.split('_')[-1] print(f'{day}: mean anomaly {v.mean():+.3f} degC')