JCSDA-internal / ufo-data

Tier-1 test files for ufo repository
1 stars 0 forks source link

Add test files for logarithm obs operator #439

Closed ReubenHill closed 1 month ago

ReubenHill commented 2 months ago

Description

Introduce files necessary for the logarithm observation operator in https://github.com/JCSDA-internal/ufo/pull/3435.

build-group=https://github.com/JCSDA-internal/ufo/pull/3435

Issue(s) addressed

N/A

Dependencies

None

Impact

None

Checklist

Code used to generate files

import netCDF4 as nc
import numpy as np

def logarithm_test_files_gen():

    missing = -3.3687952621450501176e38  # jedi's missing value

    # Descending values from 69 to 0 in the first column, 70 to 1 in the second
    # column, 71 to 2 in the third column etc
    gv = np.array([[69.0 - j + i for j in range(70)] for i in range(5)])

    ob = np.linspace(1.0, 5.0, 5)  # arbitrary
    err = np.ones(5)  # arbitrary

    nlocs = len(ob)
    nlevs = len(gv[0])
    nslevs = 1

    # Write GeoVaLs file
    filegv = nc.Dataset("logarithm_geovals.nc4", "w")
    filegv.date_time = "20200101T0000Z"
    locationgv = filegv.createDimension("nlocs", nlocs)
    nlevsgv = filegv.createDimension("nlevs", nlevs)
    nslevsgv = filegv.createDimension("nslevs", nslevs)

    latgv = filegv.createVariable("latitude", "f4", ("nlocs"))
    longv = filegv.createVariable("longitude", "f4", ("nlocs"))
    timegv = filegv.createVariable("time", "f4", ("nlocs"))
    Tgv = filegv.createVariable("air_temperature", "f4", ("nlocs", "nlevs"))
    pgv = filegv.createVariable("air_pressure", "f4", ("nlocs", "nlevs"))
    spgv = filegv.createVariable("surface_pressure", "f4", ("nlocs", "nslevs"))
    rhgv = filegv.createVariable("relative_humidity", "f4", ("nlocs", "nlevs"))

    latgv[:] = 0
    longv[:] = 0
    timegv[:] = 0
    Tgv[:] = gv[:]
    pgv[:] = gv[:] + 0.5
    spgv[:] = gv[:, -1] + 1.0
    rhgv[:] = gv[:] + 1.5

    filegv.close()

    # H(x) values (all taken at the surface)
    lnT = np.log(gv[:, -1])
    lnp = np.log(gv[:, -1] + 0.5)
    lnsp = np.log(gv[:, -1] + 1.0)
    log10T = np.log10(gv[:, -1])
    log10p = np.log10(gv[:, -1] + 0.5)
    log10sp = np.log10(gv[:, -1] + 1.0)
    log2T = np.log2(gv[:, -1])
    log2p = np.log2(gv[:, -1] + 0.5)
    log2sp = np.log2(gv[:, -1] + 1.0)
    log1p5T = np.log(gv[:, -1]) / np.log(1.5)
    log1p5p = np.log(gv[:, -1] + 0.5) / np.log(1.5)
    log1p5sp = np.log(gv[:, -1] + 1.0) / np.log(1.5)
    rh = gv[:, -1] + 1.5  # identity operator which takes the surface value

    # replace missing H(x) values with jedi's missing value
    lnT[lnT == -np.inf] = missing
    lnp[lnp == -np.inf] = missing
    lnsp[lnsp == -np.inf] = missing
    log10T[log10T == -np.inf] = missing
    log10p[log10p == -np.inf] = missing
    log10sp[log10sp == -np.inf] = missing
    log2T[log2T == -np.inf] = missing
    log2p[log2p == -np.inf] = missing
    log2sp[log2sp == -np.inf] = missing
    log1p5T[log1p5T == -np.inf] = missing
    log1p5p[log1p5p == -np.inf] = missing
    log1p5sp[log1p5sp == -np.inf] = missing
    # not necessary for rh

    # Write observation file
    fileobs = nc.Dataset("logarithm_obs.nc4", "w")
    fileobs._ioda_layout = "ObsGroup"
    fileobs._ioda_layout_version = 0
    fileobs.date_time = "20200101T0000Z"
    Location = fileobs.createDimension("Location", nlocs)

    loc = fileobs.createVariable("Location", "f4", ("Location"))
    lat = fileobs.createVariable("MetaData/latitude", "f4", ("Location"))
    lon = fileobs.createVariable("MetaData/longitude", "f4", ("Location"))
    datetime = fileobs.createVariable("MetaData/dateTime", "i8", ("Location"))
    datetime.units = "seconds since 1970-01-01T00:00:00Z"

    loc[:] = 0
    lat[:] = 0
    lon[:] = 0
    datetime[:] = 1577836800

    # Variables
    OT = fileobs.createVariable("ObsValue/airTemperature", "f4", ("Location"))
    OT[:] = ob[:]
    OP = fileobs.createVariable("ObsValue/pressure", "f4", ("Location"))
    OP[:] = ob[:] + 0.5
    OSP = fileobs.createVariable("ObsValue/surfacePressure", "f4", ("Location"))
    OSP[:] = ob[:] + 1.0
    # Errors
    OTE = fileobs.createVariable("ObsError/airTemperature", "f4", ("Location"))
    OTE[:] = err[:]
    OPE = fileobs.createVariable("ObsError/pressure", "f4", ("Location"))
    OPE[:] = err[:]
    OSPE = fileobs.createVariable("ObsError/surfacePressure", "f4", ("Location"))
    OSPE[:] = err[:]
    ORH = fileobs.createVariable("ObsValue/relativeHumidity", "f4", ("Location"))
    ORH[:] = rh[:]
    ORHE = fileobs.createVariable("ObsError/relativeHumidity", "f4", ("Location"))
    ORHE[:] = err[:]
    # Reference H(x) values
    TRLnT = fileobs.createVariable("TestReferenceLn/airTemperature", "f4", ("Location"))
    TRLnT[:] = lnT[:]
    TRLnP = fileobs.createVariable("TestReferenceLn/pressure", "f4", ("Location"))
    TRLnP[:] = lnp[:]
    TRLnSP = fileobs.createVariable(
        "TestReferenceLn/surfacePressure", "f4", ("Location")
    )
    TRLnSP[:] = lnsp[:]
    TR10T = fileobs.createVariable(
        "TestReferenceLog10/airTemperature", "f4", ("Location")
    )
    TR10T[:] = log10T[:]
    TR10P = fileobs.createVariable("TestReferenceLog10/pressure", "f4", ("Location"))
    TR10P[:] = log10p[:]
    TR10SP = fileobs.createVariable(
        "TestReferenceLog10/surfacePressure", "f4", ("Location")
    )
    TR10SP[:] = log10sp[:]
    TR2T = fileobs.createVariable(
        "TestReferenceLog2/airTemperature", "f4", ("Location")
    )
    TR2T[:] = log2T[:]
    TR2P = fileobs.createVariable("TestReferenceLog2/pressure", "f4", ("Location"))
    TR2P[:] = log2p[:]
    TR2SP = fileobs.createVariable(
        "TestReferenceLog2/surfacePressure", "f4", ("Location")
    )
    TR2SP[:] = log2sp[:]
    TR1p5T = fileobs.createVariable(
        "TestReferenceLog1p5/airTemperature", "f4", ("Location")
    )
    TR1p5T[:] = log1p5T[:]
    TR1p5P = fileobs.createVariable("TestReferenceLog1p5/pressure", "f4", ("Location"))
    TR1p5P[:] = log1p5p[:]
    TR1p5SP = fileobs.createVariable(
        "TestReferenceLog1p5/surfacePressure", "f4", ("Location")
    )
    TR1p5SP[:] = log1p5sp[:]
    # For composite, use identity for RH and log10 for T, P, SP
    TRComp = fileobs.createVariable(
        "TestReferenceComposite/relativeHumidity", "f4", ("Location")
    )
    TRComp[:] = rh[:]
    TRCompT = fileobs.createVariable(
        "TestReferenceComposite/airTemperature", "f4", ("Location")
    )
    TRCompT[:] = log10T[:]
    TRCompP = fileobs.createVariable(
        "TestReferenceComposite/pressure", "f4", ("Location")
    )
    TRCompP[:] = log10p[:]
    TRCompSP = fileobs.createVariable(
        "TestReferenceComposite/surfacePressure", "f4", ("Location")
    )
    TRCompSP[:] = log10sp[:]

    fileobs.close()

if __name__ == "__main__":
    logarithm_test_files_gen()
ReubenHill commented 2 months ago

Thank you for adding these files. Seeing the generation code is also very helpful!

I wonder if such code ought to be archived somewhere, for example in this repository.