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 [1]:

import tempfile
from pathlib import Path

import numpy as np

from earthlens import EarthLens
from earthlens.aggregate import AggregationConfig

import tempfile from pathlib import Path import numpy as np from earthlens import EarthLens from earthlens.aggregate import AggregationConfig

Request + download¶

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

In [2]:

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)')

2026-06-25 19:46:20.962 | INFO     | earthlens.erddap.backend:_fetch_grid:403 - ERDDAP griddap NOAA_DHW: GET https://coastwatch.pfeg.noaa.gov/erddap/griddap/NOAA_DHW.nc?CRW_SSTANOMALY[(2023-06-01T00:00:00Z):1:(2023-06-03T00:00:00Z)][(0.0):1:(5.0)][(150.0):1:(155.0)]

2026-06-25 19:46:51 | INFO | pyramids.base.config | Logging is configured.

CRW_SSTANOMALY_1D_20230601.tif  (12 KB)
CRW_SSTANOMALY_1D_20230602.tif  (12 KB)
CRW_SSTANOMALY_1D_20230603.tif  (11 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 [3]:

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} %')

EPSG:        4326
array shape: (101, 101)  (10,201 pixels)
valid px:    10,200
min anomaly: -0.37 degC
mean anomaly:+0.45 degC
max anomaly: +0.94 degC
pct > +1 degC (warm stress): 0.0 %

Daily area-mean anomaly across the window¶

In [4]:

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')

20230601: mean anomaly +0.445 degC
20230602: mean anomaly +0.511 degC
20230603: mean anomaly +0.589 degC