Missing "v"-component of wind from NAM GRIB files.

[bbonenfant]

I commented on #45 a few weeks ago, but I decided to create a new ticket since that ticket is closed and most likely not being tracked. For reference, I was the original reporter of that Issue. Here is the comment from that ticket:

I'm just now revisiting this, and I'm still experiencing the same bug where the u winds are present but the v winds are not. This is using cfgrib version 0.9.7.3. Below is the output from opening a NAM awphys** GRIB2 file which can be found within this database. I've made sure to check all the xarray.Dataset objects output by the cfgrib.open_datasets method and none contain v winds.

<xarray.Dataset>
Dimensions:        (isobaricInhPa: 39, x: 614, y: 428)
Coordinates:
    time           datetime64[ns] 2019-10-30
    step           timedelta64[ns] 00:00:00
  * isobaricInhPa  (isobaricInhPa) int64 1000 975 950 925 900 ... 125 100 75 50
    latitude       (y, x) float64 12.19 12.22 12.25 12.28 ... 57.39 57.36 57.33
    longitude      (y, x) float64 226.5 226.6 226.8 226.9 ... 310.3 310.4 310.6
    valid_time     datetime64[ns] 2019-10-30
Dimensions without coordinates: x, y
Data variables:
    t              (isobaricInhPa, y, x) float32 ...
    u              (isobaricInhPa, y, x) float32 ...
    w              (isobaricInhPa, y, x) float32 ...
    gh             (isobaricInhPa, y, x) float32 ...
    r              (isobaricInhPa, y, x) float32 ...
    tke            (isobaricInhPa, y, x) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             kwbc
    GRIB_centreDescription:  US National Weather Service - NCEP
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             US National Weather Service - NCEP 

Additional information requested from your contributing guidelines:

• OS: macOS High Sierra v10.13.6
• Installed using pip with xarray==0.14.1

[alexkrullwx]

Using recent RAP/RUC GRIB2 files, cfgrib has been doing that same thing with the winds. I presume the RAP is doing something similar to the winds as the NAM. My temporary work around has been to convert the GRIB2 to a netCDF4 using NCL, This does pick out the v-component from my RAP files.

Using archived HRRR files from the Utah Archive, the HRRR GRIB2 files provide both the u and v components when cfgrib puts it into an xarray for me. So far only the RAP and NAM I am seeing this with.

[alexamici]

The problem is those GRIB files use MULTI-FIELD messages for wind variables and I didn't find a way to consistently make random access reads on MULTI-FIELD messages with ecCodes.

[alexkrullwx]

The problem is those GRIB files use MULTI-FIELD messages for wind variables and I didn't find a way to consistently make random access reads on MULTI-FIELD messages with ecCodes.

Okay. Thank you for the information and update!

[alexkrullwx]

I see there have been a few updated releases in the last couple months. No specific notes on it, but has there been any luck with multi-field messages for wind in these kind of GRIB files yet?

[xebadir]

No, same issue applies. I went looking around and couldn't find a solution. NAM/RAP files seem to have the issues with the multi-fields, not an issue for GFS/HRRR GRIB 2 files.

[gdh101]

I've been struggling with this for hours, using the python eccodes library on the RAP grib2 data, and am happy to see it's not just me. Thanks for this bug report.

[xebadir]

@gdh101 I know its not a great solution, but if you need a solution for RAP use wgrib2 and just do a direct conversion to nc, you could then put it together into vertical profiles, and for an areal subset as follows (unless you are looking for direct parameters, or want the whole domain then just modify as you need):
`import cfgrib
import xarray as xr
import re
import numpy as np

#Load wgrib2 parsed file
data = xr.open_dataset('rap_filename.grb2.nc')

#Subset to domain of interest (You will want to change this)
data_sub = data.where((data.latitude>=36.)&(data.latitude<=38.)&(data.longitude>=260)&(data.longitude<=262),drop=True)

#Strip off the names of the individual variables and turn them into a list.
fieldnames = list(data_sub.keys())

#Sort the list of variables in order so that all parameters of one type will be together.
fieldnames = sorted(fieldnames, key=lambda x: x.split("_")[0])

#Search this list for all parameters containing the partial string UGRD or whichever variable.
ugrd = [x for x in fieldnames if re.search('UGRD', x)]
vgrd = [x for x in fieldnames if re.search('VGRD', x)]
temp = [x for x in fieldnames if re.search('TMP', x)]
hght = [x for x in fieldnames if re.search('HGT', x)]
rh = [x for x in fieldnames if re.search('RH', x)]

#Now Pick Out the 2D profile variables needed (may vary model to model)
ugrd_pres = [x for x in ugrd if re.search('mb', x)][0:37]
ugrd_surf = [x for x in ugrd if re.search('10maboveground', x)]
vgrd_pres = [x for x in vgrd if re.search('mb', x)][0:37]
vgrd_surf = [x for x in vgrd if re.search('10maboveground', x)]
hght_pres = hght[0:37]
temp_pres = [x for x in temp if re.search('mb', x)][0:37]
temp_surf = [x for x in temp if re.search('2maboveground', x)]
rh_pres = [x for x in rh if re.search('mb', x)][0:37]
rh_surf = [x for x in rh if re.search('2maboveground', x)]

#PRES_surface = data_sub.PRES_surface.values#[1,1].squeeze()
pres = np.array([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])

ds1 = xr.concat([data_sub[name] for name in ugrd_pres], 'pressure')
ds1 = ds1.assign_coords(pressure=pres)
ds2 = xr.concat([data_sub[name] for name in vgrd_pres], 'pressure')
ds2 =ds2.assign_coords(pressure=pres)
ds3 = xr.concat([data_sub[name] for name in temp_pres], 'pressure')
ds3 =ds3.assign_coords(pressure=pres)
ds4 = xr.concat([data_sub[name] for name in hght_pres], 'pressure')
ds4 = ds4.assign_coords(pressure=pres)
ds5 = xr.concat([data_sub[name] for name in rh_pres], 'pressure')
ds5 = ds5.assign_coords(pressure=pres)

data = xr.merge([ds1,ds2,ds3,ds4,ds5,data_sub.PRES_surface,data_sub[ugrd_surf],data_sub[vgrd_surf],data_sub[temp_surf],data_sub[rh_surf]])
data = data.rename({'UGRD_1000mb':'u', 'VGRD_1000mb':'v', 'UGRD_10maboveground': 'u10', 'VGRD_10maboveground':'v10',
'TMP_2maboveground':'t2m','RH_2maboveground':'rh2m','TMP_1000mb':'t', 'HGT_1000mb':'h','RH_1000mb':'rh','PRES_surface':'sp'})

#If you Want to output only:
data.to_netcdf('rap_filename.nc')
`
I'm sure there is a better way to key this, but its the easiest thing I could think of. Hope this helps save you some time.

[gdh101]

Thanks @xebadir. I actually was already using wgrib2 to convert to netcdf, then scipy.io to extract the data to a numpy array. It's faster iterating through a netcdf file than a grib2 (and the parameter names are resolved better), but not if accounting for the extra time to convert to netcdf. So I was trying to upgrade my code to use eccodes and ran into the v-component issue just with the RAP model so far (hadn't tried the NAM yet).

But thanks for the code. There's always a thousand ways to solve the same problem, and I'll see if there's parts I can use.

[alexkrullwx]

@gdh101 I know its not a great solution, but if you need a solution for RAP use wgrib2 and just do a direct conversion to nc, you could then put it together into vertical profiles, and for an areal subset as follows (unless you are looking for direct parameters, or want the whole domain then just modify as you need):
`import cfgrib
import xarray as xr
import re
import numpy as np

#Load wgrib2 parsed file
data = xr.open_dataset('rap_filename.grb2.nc')

#Subset to domain of interest (You will want to change this)
data_sub = data.where((data.latitude>=36.)&(data.latitude<=38.)&(data.longitude>=260)&(data.longitude<=262),drop=True)

#Strip off the names of the individual variables and turn them into a list.
fieldnames = list(data_sub.keys())

#Sort the list of variables in order so that all parameters of one type will be together.
fieldnames = sorted(fieldnames, key=lambda x: x.split("_")[0])

#Search this list for all parameters containing the partial string UGRD or whichever variable.
ugrd = [x for x in fieldnames if re.search('UGRD', x)]
vgrd = [x for x in fieldnames if re.search('VGRD', x)]
temp = [x for x in fieldnames if re.search('TMP', x)]
hght = [x for x in fieldnames if re.search('HGT', x)]
rh = [x for x in fieldnames if re.search('RH', x)]

#Now Pick Out the 2D profile variables needed (may vary model to model)
ugrd_pres = [x for x in ugrd if re.search('mb', x)][0:37]
ugrd_surf = [x for x in ugrd if re.search('10maboveground', x)]
vgrd_pres = [x for x in vgrd if re.search('mb', x)][0:37]
vgrd_surf = [x for x in vgrd if re.search('10maboveground', x)]
hght_pres = hght[0:37]
temp_pres = [x for x in temp if re.search('mb', x)][0:37]
temp_surf = [x for x in temp if re.search('2maboveground', x)]
rh_pres = [x for x in rh if re.search('mb', x)][0:37]
rh_surf = [x for x in rh if re.search('2maboveground', x)]

#PRES_surface = data_sub.PRES_surface.values#[1,1].squeeze()
pres = np.array([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])

ds1 = xr.concat([data_sub[name] for name in ugrd_pres], 'pressure')
ds1 = ds1.assign_coords(pressure=pres)
ds2 = xr.concat([data_sub[name] for name in vgrd_pres], 'pressure')
ds2 =ds2.assign_coords(pressure=pres)
ds3 = xr.concat([data_sub[name] for name in temp_pres], 'pressure')
ds3 =ds3.assign_coords(pressure=pres)
ds4 = xr.concat([data_sub[name] for name in hght_pres], 'pressure')
ds4 = ds4.assign_coords(pressure=pres)
ds5 = xr.concat([data_sub[name] for name in rh_pres], 'pressure')
ds5 = ds5.assign_coords(pressure=pres)

data = xr.merge([ds1,ds2,ds3,ds4,ds5,data_sub.PRES_surface,data_sub[ugrd_surf],data_sub[vgrd_surf],data_sub[temp_surf],data_sub[rh_surf]])
data = data.rename({'UGRD_1000mb':'u', 'VGRD_1000mb':'v', 'UGRD_10maboveground': 'u10', 'VGRD_10maboveground':'v10',
'TMP_2maboveground':'t2m','RH_2maboveground':'rh2m','TMP_1000mb':'t', 'HGT_1000mb':'h','RH_1000mb':'rh','PRES_surface':'sp'})

#If you Want to output only:
data.to_netcdf('rap_filename.nc')
`
I'm sure there is a better way to key this, but its the easiest thing I could think of. Hope this helps save you some time.

Thanks for the additional information. For the RAP I have been using NCL to convert the grib2 files to a netCDF4. Even with compression my file sizes are increasing quite a bit, but it gets the job done when I need it.

For those that are struggling with file size using this method, I have been converting them to netCDF4 files and then taking using Python to extract from just the area that I need, and then put that into a new, smaller netCDF4 file. It's time consuming, but can greatly reduce the file size.

Thanks!

[xebadir]

No, same issue applies. I went looking around and couldn't find a solution. NAM/RAP files seem to have the issues with the multi-fields, not an issue for GFS/HRRR GRIB 2 files.

Contrary to the above GFS also has this issue - went looking deeper in the variables, and it excludes them as well on the isobaric surfaces.

[alexamici]

All, I know it took a long time but I may have found a way to access consistently all data in MULTI-FIELD messages and possibly fix the cases when cfgrib silently drops the v-components.

I needed to dive deep into ecCodes and it has not been easy, but I have just committed a tentative fix in stable/0.9.8.x branch. Whoever is still interested, please test and report back if it fixes the problem for you.

[bbonenfant]

I just tested on a NAM grib and it appears that the u and v components are now available. Thank you for the continued support of this issue. I hope this turns out to be stable and reliable!

[alexamici]

Just release version 0.9.8.2 with the fix. Closing 🤞

[xebadir]

Testing for GFS and RAP also confirms this is now resolved in 0.9.8.2. Great work on solving this.

[alexkrullwx]

Excellent. Thank you for your extra effort in doing this. It is greatly appreciated. I will test it out myself by the end of this week, but it sounds like others have had success. This will be great not to have to go through the extra step of converting to netCDF4 file.

Thank you again!