Fix support for reading rotated lat/lon grids from GRIB1 files (grid type 10)

CPKalb commented 2 years ago

Describe the Enhancement

When trying to read cosmo data, the following error is generated: ERROR : gds_to_order() -> Grid type 10 not currently supported. These files are rotated latitude and longitude grids in grib1 format. lfff03080000.gz

Time Estimate

1-2 days

Sub-Issues

Consider breaking the enhancement down into sub-issues. None needed.

Relevant Deadlines

ASAP

Funding Source

7740181 WRF Modeling Using 4D for UAE

Define the Metadata

Assignee

[x] Select engineer(s) or no engineer required
[x] Select scientist(s) or no scientist required

Labels

[x] Select component(s)
[x] Select priority
[x] Select requestor(s)

Projects and Milestone

[x] Select Repository and/or Organization level Project(s) or add alert: NEED PROJECT ASSIGNMENT label
[x] Select Milestone as the next official version or Future Versions

Define Related Issue(s)

Consider the impact to the other METplus components.

[x] METplus, MET, METdatadb, METviewer, METexpress, METcalcpy, METplotpy No impacts.

Enhancement Checklist

See the METplus Workflow for details.

[x] Complete the issue definition above, including the Time Estimate and Funding Source.
[x] Fork this repository or create a branch of develop. Branch name: feature_<Issue Number>_<Description>
[x] Complete the development and test your changes.
[x] Add/update log messages for easier debugging.
[x] Add/update unit tests.
[x] Add/update documentation.
[x] Push local changes to GitHub.
[x] Submit a pull request to merge into develop. Pull request: feature <Issue Number> <Description>
[x] Define the pull request metadata, as permissions allow. Select: Reviewer(s) and Linked issues Select: Repository level development cycle Project for the next official release Select: Milestone as the next official version
[ ] Iterate until the reviewer(s) accept and merge your changes.
[ ] Delete your fork or branch.
[ ] Close this issue.

JohnHalleyGotway commented 2 years ago

Here's the COSMO website: https://www.cosmo-model.org/content/model/documentation/grib/default.htm

Running wgrib on the sample data file provides results in:

wgrib -V -d 1 lfff03080000 
rec 1:8:date 2022032000 PS kpds5=1 kpds6=1 kpds7=0 levels=(0,0) grid=255 sfc 80hr fcst:
PS=pressure [Pa]
timerange 0 P1 80 P2 0 TimeU 1  nx 281 ny 213 GDS grid 10 num_in_ave 0 missing 0
center 78 subcenter 255 process 138 Table 2 scan: WE:SN winds(grid) 
rotated LatLon grid  lat 22.000000 to 27.300000  lon 50.500000 to 57.500000
nxny 59853  (281 x 213)  dx 0 dy 0  scan 64  mode 8
transform: south pole lat -90.000000 lon 0.000000  rot angle 0.000000
min/max data 75001.5 101421  num bits 16  BDS_Ref 75001.5  DecScale 0 BinScale -1

I sure wish we actually had a sample GRIB1 file on a rotated lat/lon grid that did not have the south pole at (lat, lon) = (-90, 0). It's also curious that dx = dy = 0.

Attempting to independently plot this data with Panoply yields the following error message:

Following the log message to turn off strict GRIB1 checking leads to this next error:

So Panpoly doesn't like it.

However running ncl_convert2nc followed by ncview does produce some reasonable results. Here's surface temperature data in ncview:
The g10_lat_0 and g10_lon_1 variables showing the lat/lon values look reasonable.

The curious thing is that while this is encoded as a rotated lat/lon grid, I don't see any rotation. Might just as well be a regular lat/lon projection.

CPKalb commented 2 years ago

Thanks. I'm not as savvy with wgrib, so didn't manage to get that much information from the file. The dx and dy in python looked to be 0.025, so it seems like maybe that isn't not specified correctly in the file? Could that cause issues getting it into MET directly as a grib file?

I wondered if it actually was rotated, as it didn't seem rotated when I read it into python. I'm going to see if I can get something that looks reasonable through plot_data_plane using python embedding

JohnHalleyGotway commented 2 years ago

@CPKalb I made some local edits to the main_v10.1 branch and was able to get plot_data_plane to make this plot using that sample file.

Some details.

@rgbullock had added hooks for partially supporting Rotated Lat/Lon grids in GRIB1, but lacking sample data, they were commented out.
The GRIB1 file either explicitly specifies the dx/dy spacing or indicates that it should be derived from the min/max lat/lon values. This file indicates that they're explicitly specified... but they're specified as dx = dy = 0! And that's why wgrib reports them as such nxny 59853 (281 x 213) dx 0 dy 0 scan 64 mode 8. When I hard-code it to derive from the min/max lat/lon, I do get dx = dy = 0.25.
I could add logic so that if dx and dy are not explicitly specified as being > 0, then derive it from the min/max lat/lon... I assume that's what python is doing.

I wish we had a much wider variety of GRIB1 files on rotated lat/lon grids with which to test!

CPKalb commented 2 years ago

I'm not sure if that is what python is doing. I just printed the lat/lon values and could see that the spacing was 0.025. Probably it would be good to have logic to get dx and dy in the case that they are 0. I agree that more data would be better.

I guess the question on my end would be should I go ahead and use python embedding to read these files for uae, or is this addition to MET close to being ready?

JohnHalleyGotway commented 2 years ago

I can add the logic for accommodating this specific grid pretty easily, but I'm very concerned about the underlying assumptions. I can't imagine that these changes are correct.

They'd basically just hard code the south pole location as (lat, lon) = (-90, 0) and the rotation angle as 0.0. That's my concern. It may work for this one grid and no others, but we can't really confirm without more sample data files.

So is it a good idea to create a "partial bugfix" for a specific project/grid even though you suspect it isn't a sufficient solution?

JohnHalleyGotway commented 2 years ago

How about it make my proposed hack of a fix on a feature branch and compile it on seneca for you to test out.

While you're doing that, I'll see if I can use cnvgrib and/or wgrib to generate additional sample GRIB1 files on rotated lat/lon grids. The trouble is that I don't see in the GRIB1 spec how a rotation would even be specified. I could also inspect the wgrib source code for clues as to if/how it's done there.

CPKalb commented 2 years ago

Yes, that will work! Thanks

JohnHalleyGotway commented 2 years ago

@CPKalb you can find the bugfix branch for testing in seneca:/d1/personal/johnhg/MET/MET_development/MET-bugfix_2118_main_v10.1_grib1_rotll/met/bin. It's compiling now.

I tested this on my laptop in 2 ways:

Plotting the 2m TMP data directly:

bin/plot_data_plane lfff03080000.grib lfff03080000.ps 'name="TMP"; level="Z2";' -v 4

Regridding that field to a global lat/lon grid and then plotting the result:

bin/regrid_data_plane lfff03080000.grib G004 lfff03080000_G004.nc -field 'name="TMP"; level="Z2";' -v 10 -method NEAREST
bin/plot_data_plane lfff03080000_G004.nc lfff03080000_G004.ps 'name="TMP_Z2"; level="(*,*)";'

Here's the resulting images.

They look reasonable to me.

dtcenter / MET