Possible data parsing bugs

[RavenV-tiff]

I love the package, I'm glad someone has finally put together a weather sounding package that is easy to use and full of so many features.

There are a few issues I've noticed with how the data that gets parsed by the program. I've only tested RAP 13 km model data, so I'm not sure if this applies to any other data.

First the height calculations don't appear to be right. Through grib messages, RAP stores heights in geopotential heights (pretty standard) which should be multiplied by gravity constant to get geopotentials and then you can use the metpy function to get heights above MSL. I'm not for sure, but sounderpy may be just taking the raw geopotential heights without converting them which leads to the labeling of heights being funky (see screenshots).

Other fields look off as well. For comparison, I parsed the same RAP file with my own script and fed that into sounderpy to see if maybe it's something with how it gets plotted. However the sharppy and sounderpy plots look almost identical which leads me to believe it's the parsing.

And last minor point. The freezing level should come from the interpolated height from the sharppy prof instead of closest index to zero (or a dedicated value from the model data). There's a mismatch with the temperature trace 0C and where it gets indicated by the arrow on the plot which becomes really noticeable if plotting low vertical resolution data.

Here's the script I used to select data from RAP to compare with so maybe there's an issue with that also I'm not catching.

import numpy as np
import csv
import sharppy
import sharppy.sharptab.profile as profile
import sharppy.sharptab.interp as interp
import sharppy.sharptab.winds as winds
import sharppy.sharptab.utils as utils 
import sharppy.sharptab.params as params
import sharppy.sharptab.thermo as thermo
import sharppy.sharptab.constants as const
import pygrib
import metpy
import metpy.calc as mpcalc
from metpy.units import units
import sounderpy as spy
from sounderpy.calc import sounding_params
import requests
###OPEN FILE###
target_lat =33.84
target_lon =-102.72
grib = pygrib.open('rap_130_20220524_0000_000.grb2')
###GET ELEVATION###
def get_elevation(lat, lon):
    base_url = "https://api.open-elevation.com/api/v1/lookup"
    params = {
        "locations": f"{lat},{lon}",
        "key": "value"  # Dummy key, not required for Open Elevation API
    }
    response = requests.get(base_url, params=params)
    data = response.json()
    if "results" in data and data["results"]:
        elevation = data["results"][0]["elevation"]
        return elevation
    else:
        return None

elevation = get_elevation(target_lat, target_lon)
###EXTRACT MODEL DATA###
def find_nearest_latlon(grib, target_lat, target_lon):
    """
    Find the nearest GRIB message to the specified latitude and longitude.
    """
    result_list = []  # List to store results for each message
    for grb in grib:
        latitudes, longitudes = grb.latlons()
        distances = ((latitudes - target_lat)**2 + (longitudes - target_lon)**2)**0.5
        min_index = distances.argmin()
        nearest_grb_value = grb.values.flat[min_index]
        nearest_grb = {
            'message_number': grb.messagenumber,
            'name': grb.name,
            'latitude': latitudes.flat[min_index],
            'longitude': longitudes.flat[min_index],
            'distance': distances.flat[min_index],
            'value': nearest_grb_value
        }
        result_list.append(nearest_grb)
    return result_list

result_list = find_nearest_latlon(grib, target_lat, target_lon)
###PARSE RAP DATA###
pres = [1000,975,950,925,900,875,850,825,800,775,750,725,700,675,650,625,600,575,550,525,500,475,450,425,400,375,350,325,300,275,250,225,200,175,150,125,100]
surface_pressure = result_list[232]['value'] / 100
index_array = [231,219,212,206,200,194,188,182,176,170,164,158,152,146,140,134,128,122,116,110,103,97,91,85,79,73,67,61,55,49,43,37,31,25,19,13,7]
values = [result_list[idx]['value'] for idx in index_array]
hght = [mpcalc.geopotential_to_height(units.Quantity(value*9.80665, 'm^2/s^2')).m for value in values]
#hght = list(map(lambda x: x - elevation, hght))
index_array = [225,220,213,207,201,195,189,183,177,171,165,159,153,147,141,135,129,123,117,111,104,98,92,86,80,74,68,62,56,50,44,38,32,26,20,14,8]
values = [result_list[idx]['value'] for idx in index_array]
tmpc = [(value*units.degK).to(units.degC).m for value in values]
index_array = [226,221,214,208,202,196,190,184,178,172,166,160,154,148,142,136,130,124,118,112,105,99,93,87,81,75,69,63,57,51,45,39,33,27,21,15,9]
values = [result_list[idx]['value'] for idx in index_array]
dwpc = [(mpcalc.dewpoint_from_relative_humidity(units.Quantity(t, 'degC'), units.Quantity(v, 'percent')).m, v)[0] for t, v in zip(tmpc, values)]
index_array = [228,223,216,210,204,198,192,186,180,174,168,162,156,150,144,138,132,126,120,114,107,101,95,89,83,77,71,65,59,53,47,41,35,29,23,17,11]
values = [result_list[idx]['value'] for idx in index_array]
index_array2 = [229,224,217,211,205,199,193,187,181,175,169,163,157,151,145,139,133,127,121,115,108,102,96,90,84,78,72,66,60,54,48,42,36,30,24,18,12]
values2 = [result_list[idx]['value'] for idx in index_array2]
wspd = [(mpcalc.wind_speed(units.Quantity(u, 'm/s'), units.Quantity(v, 'm/s')).to(units.knot).m) for u, v in zip(values, values2)]
wdir = [(mpcalc.wind_direction(units.Quantity(u, 'm/s'), units.Quantity(v, 'm/s')).m.item()) for u, v in zip(values, values2)]

###REMOVE SUBSUFACE###
hght = [elevation] + hght[next(i for i, p in enumerate(pres) if p < surface_pressure):]
tmpc = [(result_list[239]['value']*(units.degK)).to(units.degC).m] + [tmp for tmp, p in zip(tmpc, pres) if p < surface_pressure]
dwpc = [(result_list[242]['value']*(units.degK)).to(units.degC).m] + [tmp for tmp, p in zip(dwpc, pres) if p < surface_pressure]
wspd = [mpcalc.wind_speed(units.Quantity(result_list[247]['value'], 'm/s'), units.Quantity(result_list[248]['value'], 'm/s')).to(units.knot).m] + [tmp for tmp, p in zip(wspd, pres) if p < surface_pressure]
wdir = [mpcalc.wind_direction(units.Quantity(result_list[247]['value'], 'm/s'), units.Quantity(result_list[248]['value'], 'm/s')).m.item()] + [tmp for tmp, p in zip(wdir, pres) if p < surface_pressure]
pres = [surface_pressure] + [level for level in pres if level < surface_pressure]
###SHARPPY PROFILE###
prof = profile.create_profile(profile='default', pres=pres, hght=hght, tmpc=tmpc, \
                                    dwpc=dwpc, wspd=wspd, wdir=wdir, missing=-9999, strictQC=True)

###SPY DICTIONARY CREATION###
data2 = {}
data2['T'] = tmpc*units.degC
data2['Td'] = dwpc*units.degC
data2['rh'] = (mpcalc.relative_humidity_from_dewpoint(data2['T'], dwpc*units.degC)*100)*units.percent
index_array = [228,223,216,210,204,198,192,186,180,174,168,162,156,150,144,138,132,126,120,114,107,101,95,89,83,77,71,65,59,53,47,41,35,29,23,17,11]
values = [result_list[idx]['value'] for idx in index_array]
data2['u'] = ([result_list[247]['value']] + [tmp for tmp, p in zip(values, pres) if p < surface_pressure])*units('m/s').to(units.knot)
index_array = [229,224,217,211,205,199,193,187,181,175,169,163,157,151,145,139,133,127,121,115,108,102,96,90,84,78,72,66,60,54,48,42,36,30,24,18,12]
values = [result_list[idx]['value'] for idx in index_array]
data2['v'] = ([result_list[248]['value']] + [tmp for tmp, p in zip(values, pres) if p < surface_pressure])*units('m/s').to(units.knot)
data2['z'] = list(map(lambda x: x - elevation, hght))*units('meters')
data2['p'] = pres*units('hectopascal')
data2['zAGL'] = hght*units('meters')
data2['site_info'] = data['site_info']

spy.build_sounding(data2, save=True, filename='test')

Screenshots Original Sounderpy plot when pulling data from NCEI Sounderpy plot when fed my parsed RAP file Sharppy plot when exporting my parsed RAP file Cameron Nixon's plot, which looks similar to my parsed plot, though I believe he's modified some of the original model values?

[kylejgillett]

Hey @RavenV-tiff!

Thanks a lot for opening this issue. Upon further investigation, it turns out that there was an issue with the way I was creating the pressure array after finding the true surface. I fixed the weird height data, and in the process opted to change the way SounderPy parsed out specific lat/lon information. Instead of picking a point, SounderPy will use a 'box average' approach, which can be user-specified. The default is 0.1 degrees.

I also changed the frz point calculation to an interpolation scheme which is far more accurate. Nice catch.

A bug-fix release is on the way soon.

[RavenV-tiff]

Very nice! Look forward to the next release :)