getMH1SST_Sci.py
Science quality MODIS-Aqua SST retrieval and qual_sst masking.
"""
Retrieve, mask, and publish science-quality MODIS-Aqua SST products.
This script queries the NASA OceanColor file search API for MODIS-Aqua
Level-3 mapped science-quality sea surface temperature (SST) products,
compares the returned filenames against files already staged in the production
Google Cloud Storage (GCS) bucket, downloads only missing files to the Cloud
Run container, publishes the raw SST files, creates masked SST files,
publishes the masked files, and removes local temporary files.
For each downloaded SST file, the script creates a masked variant in which
pixels with qual_sst < 0 are written as fill values. Pixels with qual_sst >= 0
are retained in the output variable sstMasked.
Periods
-------
- DAY -> 1day
- 8D -> 8day
- MO -> mday
NASA query settings
-------------------
- sensor_id: 7 (MODIS-Aqua)
- dtid: 1049
- resolution_id: 4km
- stream: science quality / standard delayed
Output
------
Raw SST files are published under the MH1 SST path family. Masked SST files
are published under the corresponding MH1 sstMask path family. Science-quality
paths do not use the _NRT suffix.
Runtime requirements
--------------------
- ROYLIB_CONFIG points to the runtime config.yml file.
- Earthdata .netrc and URS cookies are available under /tmp.
- The Cloud Run service account can read and write the configured GCS bucket.
- netCDF4 and numpy are available in the container.
Usage
-----
python getMH1SST_Sci.py
"""
if __name__ == "__main__":
from datetime import datetime
import os
import subprocess
from pathlib import Path
from netCDF4 import Dataset, date2num
import numpy as np
import numpy.ma as ma
from roylib import CFG, list_bucket_content, send_to_servers
# --- Config Validation ---
for _k in ("ERDPROD_BUCKET", "HOME_DIR", "NASA_FILE_SEARCH_URL", "NASA_GETFILE_URL"):
if not CFG.get(_k):
raise KeyError(f"CFG['{_k}'] is required but not set")
# --- Constants ---
# dtid=1049: MODIS-Aqua SST Standard Science Quality L3 mapped
DTID = "1049"
periods = {
"DAY": {"dir": "1day", "flag": "1"},
"8D": {"dir": "8day", "flag": "8"},
"MO": {"dir": "mday", "flag": "m"},
}
now = datetime.now()
start_date = f"{now.year}-01-01"
end_date = now.strftime("%Y-%m-%d")
bucket_name = CFG.get("ERDPROD_BUCKET")
home_dir = Path(CFG.get("HOME_DIR", "/tmp"))
cookie_path = home_dir / ".urs_cookies"
download_dir = home_dir / "data"
resources_dir = home_dir / "resources"
download_dir.mkdir(parents=True, exist_ok=True)
resources_dir.mkdir(parents=True, exist_ok=True)
# -----------------------------------------------------------------------
def _make_masked_nc(src_path: Path, period_flag: str) -> Path | None:
"""
Create a masked SST NetCDF file from a raw MODIS-Aqua SST file.
The raw SST variable is read from src_path. Pixels where qual_sst < 0
are masked and written as fill values. Pixels where qual_sst >= 0 are
retained in the output variable sstMasked.
Parameters
----------
src_path : pathlib.Path
Path to the raw SST NetCDF file.
period_flag : str
Product period flag used by send_to_servers(); expected values are
"1", "8", or "m".
Returns
-------
pathlib.Path | None
Path to the masked NetCDF file if creation succeeds; otherwise None.
"""
try:
with Dataset(str(src_path), "r") as ds:
sst_raw = ds.variables["sst"][:, :]
qual = ds.variables["qual_sst"][:, :]
global_attrs = {a: getattr(ds, a) for a in ds.ncattrs()}
lat_vals = ds.variables["lat"][:]
lon_vals = ds.variables["lon"][:]
except Exception as exc:
print(f"[WARN] Could not read {src_path.name} for masking: {exc}")
return None
sst_masked = ma.array(sst_raw, mask=(qual < 0), fill_value=-999.0)
stem = src_path.stem
masked_name = stem.replace(".sst.", ".sstMasked.") + ".nc"
masked_path = download_dir / masked_name
# Derive centre time from filename date stamp (AQUA_MODIS.YYYYMMDD...)
try:
fname = src_path.name
y, m, d = int(fname[11:15]), int(fname[15:17]), int(fname[17:19])
if period_flag == "1":
centre = datetime(y, m, d, 12, 0, 0)
time_val = date2num(centre, units="seconds since 1970-01-01")
elif period_flag == "8":
y2, m2, d2 = int(fname[20:24]), int(fname[24:26]), int(fname[26:28])
t1 = date2num(datetime(y, m, d), units="seconds since 1970-01-01")
t2 = date2num(datetime(y2, m2, d2), units="seconds since 1970-01-01")
time_val = (t1 + t2) / 2.0
else: # monthly
centre = datetime(y, m, 16, 0, 0)
time_val = date2num(centre, units="seconds since 1970-01-01")
except Exception:
time_val = date2num(now, units="seconds since 1970-01-01")
try:
with Dataset(str(masked_path), "w", format="NETCDF4_CLASSIC") as out:
out.createDimension("lat", len(lat_vals))
out.createDimension("lon", len(lon_vals))
out.createDimension("time", 1)
for attr, val in global_attrs.items():
try:
setattr(out, attr, val)
except Exception:
pass
out.source_data = src_path.name
out.creation_date = str(now)
out.mask_applied = "qual_sst >= 0"
t_var = out.createVariable("time", "f8", ("time",))
t_var.units = "seconds since 1970-01-01"
t_var[0] = time_val
lat_var = out.createVariable("lat", "f4", ("lat",))
lat_var.units = "degrees_north"
lat_var[:] = lat_vals
lon_var = out.createVariable("lon", "f4", ("lon",))
lon_var.units = "degrees_east"
lon_var[:] = lon_vals
sst_var = out.createVariable(
"sstMasked", "f4", ("time", "lat", "lon"),
fill_value=-999.0, zlib=True, complevel=2
)
sst_var.long_name = "Sea Surface Temperature (qual_sst >= 0 mask applied)"
sst_var.units = "degree_C"
sst_var[0, :, :] = sst_masked[:, :]
except Exception as exc:
print(f"[WARN] Could not write masked file {masked_name}: {exc}")
return None
return masked_path
# -----------------------------------------------------------------------
for nasa_per, per_info in periods.items():
# Science Quality omits the _NRT suffix in the GCS path
stage_dir_raw = f"satellite/MH1/sst/{per_info['dir']}"
raw_staged = list_bucket_content(bucket_name, stage_dir_raw)
staged_flist = [os.path.basename(f).strip() for f in raw_staged] if raw_staged else []
print(f"\n[START] Science Quality SST {nasa_per} -> {stage_dir_raw}")
# --- Query NASA file_search API ---
query_params = (
f"results_as_file=1&sensor_id=7&dtid={DTID}&subType=1&resolution_id=4km"
f"&sdate={start_date} 00:00:00&edate={end_date} 23:59:59"
f"&prod_id=sst&period={nasa_per}"
)
query_url = f'{CFG["NASA_FILE_SEARCH_URL"]} --post-data="{query_params}"'
query_out = resources_dir / f"sst_sci_{nasa_per}_list.txt"
subprocess.run(
f"wget -4 --no-check-certificate -O {query_out} {query_url}",
shell=True, capture_output=True
)
if not query_out.exists():
print(f"[WARN] Query output not found for {nasa_per}, skipping.")
continue
# --- Download & process missing files ---
with open(query_out, "r") as f:
for line in f:
fname = line.strip()
if not fname or "No Results Found" in fname or "<html" in fname.lower():
continue
if fname in staged_flist:
continue
out_path = download_dir / fname
print(f"[DOWNLOAD] {fname}")
wget_cmd = (
f'wget -4 --netrc --auth-no-challenge=on --keep-session-cookies '
f'--load-cookies {cookie_path} --save-cookies {cookie_path} '
f'--content-disposition --no-check-certificate '
f'-O "{out_path}" {CFG["NASA_GETFILE_URL"]}/{fname}'
)
result = subprocess.run(wget_cmd, shell=True, capture_output=True)
if result.returncode != 0:
print(f"[WARN] wget failed for {fname}, skipping.")
continue
# -- Upload raw SST --
try:
send_to_servers(str(out_path), "MH1/sst", per_info["flag"])
except Exception as exc:
print(f"[ERROR] Raw upload failed for {fname}: {exc}")
# -- Build and upload masked SST --
masked_path = _make_masked_nc(out_path, per_info["flag"])
if masked_path and masked_path.exists():
try:
send_to_servers(str(masked_path), "MH1/sstMask", per_info["flag"])
except Exception as exc:
print(f"[ERROR] Masked upload failed for {masked_path.name}: {exc}")
finally:
masked_path.unlink(missing_ok=True)
# -- Cleanup --
out_path.unlink(missing_ok=True)
print("\n[COMPLETE] Science Quality SST Suite finished.")