tjturnage / cloud-radar-server

Radar server in the cloud for NOAA computing project
MIT License
2 stars 2 forks source link

Verify placefile shifting is working correctly #24

Closed lcarlaw closed 2 months ago

lcarlaw commented 5 months ago

The 5/20/2013 NSE placefiles don't look correct.

Verify the NSE code and subsequent placefile output is working as expected. Initial inspection of the data suggests that the incoming RAP files may be formatted in a different manner resulting in output that's incorrect. These files were downloaded from the NCEI archives and are pressure-level grids instead of hybrid-sigma, although logic in the reading functions does account for this difference.

shifting_error

lcarlaw commented 5 months ago

PR https://github.com/tjturnage/cloud-radar-server/pull/26 addresses (hopefully) the jagged edges which were caused by geojsoncontour skipping among individual paths. Shifting appears to be working as expected, at least for MLCAPE which was tested.

TODO: look at the other stp components (mllcl, esrh, ebwd, and mlcin) to verify calculations are working as expected.

lcarlaw commented 3 months ago

I believe this has been verified as working. There appears to be something up with the actual model data used for this event which is causing errors with the effective inflow base and top.

lcarlaw commented 2 months ago

This issue cropped up again in recent testing.

Finally tracked the issue down to a change in how more recent versions of numba implicitly handle try-except statements, which was causing the wrong part of the wobf function in sharptab/thermo.py to execute. A minimal reproduction of the issue is shown below. In short:

import numpy as np
from numba import njit 

@njit
def wobf(t):
    t = t - 20
    try:
        # If t is a scalar
        if t <= 0:
            print('t <= 0 block')
            npol = 1. + t * (-8.841660499999999e-3 + t * ( 1.4714143e-4 + t * (-9.671989000000001e-7 + t * (-3.2607217e-8 + t * (-3.8598073e-10)))))
            npol = 15.13 / (np.power(npol,4))
            return npol
        else:
            print('t > 0 block')
            ppol = t * (4.9618922e-07 + t * (-6.1059365e-09 + t * (3.9401551e-11 + t * (-1.2588129e-13 + t * (1.6688280e-16)))))
            ppol = 1 + t * (3.6182989e-03 + t * (-1.3603273e-05 + ppol))
            ppol = (29.93 / np.power(ppol,4)) + (0.96 * t) - 14.8
            return ppol

    except:
        print('--- array block ---')
        # If t is an array
        npol = 1. + t * (-8.841660499999999e-3 + t * ( 1.4714143e-4 + t * (-9.671989000000001e-7 + t * (-3.2607217e-8 + t * (-3.8598073e-10)))))
        npol = 15.13 / (np.power(npol,4))
        ppol = t * (4.9618922e-07 + t * (-6.1059365e-09 + t * (3.9401551e-11 + t * (-1.2588129e-13 + t * (1.6688280e-16)))))
        ppol = 1 + t * (3.6182989e-03 + t * (-1.3603273e-05 + ppol))
        ppol = (29.93 / np.power(ppol,4)) + (0.96 * t) - 14.8
        correction = np.zeros(t.shape, dtype=np.float64)
        correction[t <= 0] = npol[t <= 0]
        correction[t > 0] = ppol[t > 0]
        return correction

t = np.array([-1, 0, 1, 2, 3], dtype='float64')
print(f"wobf function output:\n"
    f"{wobf(t)}")

print(
    f"The correct values should be:\n"
    f"[6.103814   6.41105332 6.73092552 7.06360765 7.40923593]")

This commit attempts to work around this issue. Data for the wobf function is always evaluated as an array, with a separate wobf function created to handle scalar instances.

An example of the errors in the effective inflow layer (specifically the mu-parcel's top pressure) are below:

Error using original code (valid: 2021-12-11 04z) esrh_error

Following corrections: esrh

Matches much more favorably with the SPC mesoanalysis archive: spc-meso