diff --git a/scripts/nldas/process_nldasv2_noah.py b/scripts/nldas/process_nldasv2_noah.py index 08d34df0cb..191cc7f93b 100644 --- a/scripts/nldas/process_nldasv2_noah.py +++ b/scripts/nldas/process_nldasv2_noah.py @@ -8,31 +8,13 @@ import subprocess import sys -import pygrib from pyiem.iemre import hourly_offset from pyiem.util import ncopen, utc -SOIL_LEVELS = { - "0-10 cm down": 0, - "10-40 cm down": 1, - "40-100 cm down": 2, - "100-200 cm down": 3, -} -DEPTHS = { - "0-10 cm down": 100.0, - "10-40 cm down": 300.0, - "40-100 cm down": 600.0, - "100-200 cm down": 1000.0, -} -VAL2NC = { - "TSOIL": "soilt", - "SOILM": "soilm", -} - def process(valid): """Run for the given UTC timestamp.""" - fn = f"nldas.{valid:%Y%m%d%H}" + fn = f"nldas.{valid:%Y%m%d%H}.nc" cmd = [ "wget", "-q", @@ -44,40 +26,37 @@ def process(valid): "-O", fn, "https://hydro1.gesdisc.eosdis.nasa.gov/data/NLDAS/" - f"NLDAS_NOAH0125_H.002/2023/{valid:%03j}/NLDAS_NOAH0125_H." - f"A{valid:%Y%m%d}.{valid:%H}00.002.grb", + f"NLDAS_NOAH0125_H.2.0/{valid:%Y}/{valid:%03j}/NLDAS_NOAH0125_H." + f"A{valid:%Y%m%d}.{valid:%H}00.020.nc", ] if not os.path.isfile(fn): subprocess.call(cmd) - # pygrib can't deal with this grib1 file, so we have to hack away here - with subprocess.Popen(["wgrib", fn], stdout=subprocess.PIPE) as proc: - content = proc.stdout.read().decode("utf-8") - grbs = pygrib.open(fn) idx = hourly_offset(valid) - for line in content.split("\n"): - tokens = line.split(":") - if len(tokens) < 11: - continue - if tokens[3] not in ["TSOIL", "SOILM"]: - continue - if tokens[11] not in SOIL_LEVELS: - continue - grb = grbs[int(tokens[0])] - if tokens[3] == "SOILM": - vals = grb.values / DEPTHS[tokens[11]] - else: - vals = grb.values - with ncopen(f"/mesonet/data/nldas/{valid:%Y}_hourly.nc", "a") as nc: - nc.variables[VAL2NC[tokens[3]]][ - idx, SOIL_LEVELS[tokens[11]] - ] = vals + SOILTVARS = [ + "SoilT_0_10cm", + "SoilT_10_40cm", + "SoilT_40_100cm", + "SoilT_100_200cm", + ] + DEPTH_MM = [100.0, 300.0, 600.0, 1000.0] + with ncopen(fn) as ncin, ncopen( + f"/mesonet/data/nldas/{valid:%Y}_hourly.nc", "a" + ) as ncout: + for i, ncvar in enumerate(SOILTVARS): + ncout.variables["soilt"][idx, i] = ncin.variables[ncvar][0] + # mm over depth, we want fraction + ncinv = ncvar.replace("SoilT", "SoilM") + ncout.variables["soilm"][idx, i] = ( + ncin.variables[ncinv][0] / DEPTH_MM[i] + ) + os.unlink(fn) def main(argv): """Run for a given UTC date.""" valid = utc(int(argv[1]), int(argv[2]), int(argv[3])) - for hr in range(1): + for hr in range(24): process(valid.replace(hour=hr))