search

ufs-community / UFS_UTILS

Utilities for the NCEP models.
Other
21 stars 104 forks source link

Can we remove wgrib2 from the dependencies? #591

Closed edwardhartnett closed 2 years ago

edwardhartnett commented 2 years ago

Somehow wgrib2 crept into our dependencies.

We are trying to eliminate use of wgrib2 in NCEPLIBS. The wgrib2 library is not well-maintained, distributed, or tested, and there is no interest on the wgrib2 team in meeting our needs in these areas. So we need to not build more tools on top of wgrib2.

kgerheiser commented 2 years ago

UFS_UTILS definitely depends on the wgrib2 library. That hasn't changed recently. The alternative is to port UFS_UTILS to use g2 as its grib2 library.

GeorgeGayno-NOAA commented 2 years ago

chgres_cube uses wgrib2. It could be modified to use the 'g2' library. I don't think that would be difficult. But it would take some work and testing. And it could be done is smaller chunks.

edwardhartnett commented 2 years ago

@GeorgeGayno-NOAA for future work, can you communicate to the other chgres_cube programmers that wgrib2 should be avoided?

GeorgeGayno-NOAA commented 2 years ago

@GeorgeGayno-NOAA for future work, can you communicate to the other chgres_cube programmers that wgrib2 should be avoided?

Absolutely.

Maybe I can add a warning to our wiki page? https://github.com/NOAA-EMC/UFS_UTILS/wiki Or the PR template? https://github.com/NOAA-EMC/UFS_UTILS/blob/develop/.github/PULL_REQUEST_TEMPLATE

edwardhartnett commented 2 years ago

If we can remove the dependency on wgrib2, we should.

GeorgeGayno-NOAA commented 2 years ago

I added a warning page to the wiki: https://github.com/ufs-community/UFS_UTILS/wiki/A-Warning-about-the-WGRIB2-library

XianwuXue-NOAA commented 2 years ago

@GeorgeGayno-NOAA Based on the warning page "Pull requests that use WGRIB2 will be rejected.", do it also apply to the GEFS WCOSS2 Transition? Because GEFS uses chgres_cube.

GeorgeGayno-NOAA commented 2 years ago

@GeorgeGayno-NOAA Based on the warning page "Pull requests that use WGRIB2 will be rejected.", do it also apply to the GEFS WCOSS2 Transition? Because GEFS uses chgres_cube.

Are you adding new code to chgres_cube that uses wgrib2?

XianwuXue-NOAA commented 2 years ago

No. I didn't change code except changing "grb2_UNDEFINED" to "9.999e20" to use wgrib2/2.0.7. My changes are https://github.com/XianwuXue-NOAA/UFS_UTILS/compare/1b44601...a825b29

GeorgeGayno-NOAA commented 2 years ago

Count number of soil layers using g2 instead of wgrib2 - d327b7a. All consistency tests passed on Dell.

GeorgeGayno-NOAA commented 2 years ago

Module model_grid.F90 sets up the input grid ESMF grid object. This requires the lat/lons of the input grid centers and corner points. The wgrib2 and G2 versions give the same result for GFS lat/lon grids. But for input data on lambert conformal and rotated lat/lon grids, the computed lat/lons are different. Details to follow.

GeorgeGayno-NOAA commented 2 years ago

Ran the GRIB2 regression tests on Orion using my branch at 2425c5d. Compare the input grid lat/lons against 'develop' at 04700f9 for the 13km.conus.rap.grib2.sh test. The input RAP data for that test is on a rotated lat/lon grid. Currently, in 'develop', the lat/lons are not computed using wgrib2; Rather, they are read in from file - latlon_grid3.32739.nc. The branch computes lat/lons using IP library routine gdswzd.

Output of the grid cell center lat/lons of the four corner points and the mid-point of the grid.

'develop':

 wgrib2 lat/lon 11   -10.59060      -139.0858
 wgrib2 lat/lon ni/11   -10.59060      -72.91415
 wgrib2 lat/lon 11/nj    46.59194       125.3390
 wgrib2 lat/lon ni/nj    46.59194       22.66102
 wgrib2 lat/lon mid      53.93890      -106.2070

Branch:

after gdswzd lat/lon 11  -10.5910000000000        -139.08600
after gdswzd lat/lon ni/11  -10.5906866621875        -72.91433
after gdswzd lat/lon 11/nj   46.5915517972299        125.338962048687
after gdswzd lat/lon ni/nj   46.5920000000000        22.6610000000005
after gdswzd lat/lon mid     53.9388096007483        -106.20757

Output of the grid cell corners of the four corner points and the mid-point of the grid.

'develop':

wgrib2 lat/lon corner 11      -10.66110      -139.0706
wgrib2 lat/lon corner ni/11   -10.66110      -72.92938
wgrib2 lat/lon corner 11/nj    46.55700       125.2473
wgrib2 lat/lon corner ni/nj    46.55700       22.75274
wgrib2 lat/lon corner mid      53.87777      -106.3100

Branch:

after gdswzd lat/lon corner 11  -10.6615060520909       -139.0707
after gdswzd lat/lon corner ni/11  -10.6611928420218   -72.92955
after gdswzd lat/lon corner 11/nj   46.5566223380035    125.24724795
after gdswzd lat/lon corner ni/nj   46.5570699662187     22.752714383
after gdswzd lat/lon corner mid     53.8776705307114   -106.31061

The input RAP data has an approximate resolution of 12km. So the differences between the gdswzd calculation and the values in the latlon_grid3.32739.nc. file are small.

GeorgeGayno-NOAA commented 2 years ago

Here, the lat/lons computed by 'develop' and the branch (at 2425c5d) for the 13km.conus.nam.grib2 test are compared. The input data is NAM GRIB2 data on a lambert conformal grid. The branch computes the lat/lons using IP routine gdswzd. 'develop' uses wgrib2 to compute the grid cell center lat/lons. But for the grid cell corners, 'develop' uses routine get_cell_corners.

Output of the grid cell center lat/lons of the four corner points and the mid-point of the grid.

'develop'

wgrib2 lat/lon 11       12.19000       226.5410
wgrib2 lat/lon ni/11    14.34209       294.8745
wgrib2 lat/lon 11/nj    54.56534       207.1214
wgrib2 lat/lon ni/nj    57.32843       310.5840
wgrib2 lat/lon mid      40.56840       259.3559

branch:

after gdswzd lat/lon 11      12.1900000000000        226.541000000000
after gdswzd lat/lon ni/11   14.3420196402024        294.874751547676
after gdswzd lat/lon 11/nj   54.5655092817601        207.121247173388
after gdswzd lat/lon ni/nj   57.3285185322474        310.584563462297
after gdswzd lat/lon mid     40.5685090018817        259.356124993238

For a 13 km grid, the grid cell center lat/lons match quite well.

Output of the grid cell corners of the four corner points and the mid-point of the grid.

'develop'

wgrib2 lat/lon corner 11       12.1351774975895        226.597069691178
wgrib2 lat/lon corner ni/11    14.2872627161717        294.817884713184
wgrib2 lat/lon corner 11/nj    54.6201214681951        207.216058832940
wgrib2 lat/lon corner ni/nj    57.3832094588327        310.482317962108
wgrib2 lat/lon corner mid      40.5135624856651        259.428029461158

branch:

after gdswzd lat/lon corner 11      12.1236721352827   226.503794008511
after gdswzd lat/lon corner ni/11   14.2776458634527   294.916856390839
after gdswzd lat/lon corner 11/nj   54.5887737437851   207.013987601792
after gdswzd lat/lon corner ni/nj   57.3563259130907   310.691578562183
after gdswzd lat/lon corner mid     40.5136305429862   259.289718875438

Here, the differences seem large. And note that the corner longitude computed by 'develop' are located to the east of the grid cell center longitude. Example: at point (1,1) 226.597 is east of 226.541. Something does not seem correct.

A similar result is seen in the regression tests that use the HRRR GRIB2 data as input, which is also on a lambert conformal grid.

kgerheiser commented 2 years ago

@GeorgeGayno-NOAA what you think the issue is? Is it g2 or gdswzd?

GeorgeGayno-NOAA commented 2 years ago

@GeorgeGayno-NOAA what you think the issue is? Is it g2 or gdswzd?

We don't know yet. I suspect the problem is in wgrib2 - not g2 or gdswzd.

GeorgeGayno-NOAA commented 2 years ago

Some further analysis of https://github.com/ufs-community/UFS_UTILS/issues/591#issuecomment-1023656002 and https://github.com/ufs-community/UFS_UTILS/issues/591#issuecomment-1024381241

My branch at 474401f and 'develop' at 26cd024 were modified to output the grid cell center and corner lat/lons in netcdf format.

The files were placed on Hera in /scratch1/NCEPDEV/da/George.Gayno/ufs_utils.git/UFS_UTILS.upstream/latlon_check.

The latlon.gdswzd.nc files were produced by my branch and the latlon.wgrib2.nc files were produced by 'develop' using wgrib2. Note, lat/lons for the 13km RAP grid are not created from wgrib2; they are read in from file. That file is latlon_grid3.32769.nc. The ncdiff utility was used to create difference files (diff.nc). And the nccmp utility was used to print a summary of the differences (diff.log).

The files are in these sub-directories:

More analysis will follow.

GeorgeGayno-NOAA commented 2 years ago

Here are the differences between my branch and 'develop' for the 13km RAP grid.

Variable        Group  Count         Sum      AbsSum          Min         Max      Range         Mean      StdDev
gridlat         /     790759    -120.266     162.903 -0.000530243 0.000289917 0.00082016 -0.000152089 0.000204953
gridlon         /     792752       49.76     771.789     -79.9549         360    439.955  6.27686e-05    0.414187
gridlat_corners /     792644    -120.655     164.868 -0.000541687  0.00219727 0.00273895 -0.000152219 0.000206831
gridlon_corners /     794497    -229.414     323.772    -0.181824    0.178787   0.360611 -0.000288754 0.000780303

The mean differences are very small - much smaller than the resolution of the grid. Viewing the files in ncview show the large longitude differences (-79.9 and 360 degrees) were located at the north pole (where the longitude lines converge) and the dateline respectively (where one point was -180 and the other was 180.).

Since we know the latlon_grid3.32769.nc file is correct, we can conclude the branch values are correct as well.

GeorgeGayno-NOAA commented 2 years ago

Here are the differences between the branch 474401f and 'develop' for the 13km NAM grid (lambert conformal)

Variable        Group  Count      Sum  AbsSum          Min         Max       Range        Mean      StdDev
gridlat         /     261655  22.6143 23.3464 -6.58035e-05 0.000183105 0.000248909 8.64278e-05 5.58303e-05
gridlon         /     254008  48.8392 50.8008 -0.000137329 0.000549316 0.000686646 0.000192274 0.000141188
gridlat_corners /     263804  513.784 2193.29   -0.0313454   0.0270767   0.0584221   0.0019476  0.00970519
gridlon_corners /     263835 -36764.1 36866.6    -0.212463    0.209259    0.421722   -0.139345   0.0287544

The grid cell centers are very close. The differences are much smaller than the resolution of the grid. The grid cell corners are off - especially the longitudes. The branch computes the centers using wgrib2 and the corners using routine get_cell_corners.

A similar result is seen with the 3km HRRR grid. The corner longitude difference is similar to the resolution of the grid:

Variable        Group   Count      Sum  AbsSum          Min          Max       Range         Mean      StdDev
gridlat         /     1885229   98.535 149.762 -0.000123978  0.000217438 0.000341415  5.22669e-05  7.8939e-05
gridlon         /     1905141 -871.205 871.205 -0.000732422 -0.000244141 0.000488281 -0.000457291 0.000110639
gridlat_corners /     1907294  830.204 4353.38  -0.00643921    0.0065918    0.013031  0.000435279   0.0026419
gridlon_corners /     1908000 -65904.6 65961.1   -0.0474243    0.0458069   0.0932312   -0.0345412  0.00527456
GeorgeGayno-NOAA commented 2 years ago

The second - and hopefully - final PR was created for this issue - #641. This PR replaces the wgrib2 computation of input grid latitude/longitude with a computation from IPLIB routine gdswzd. So, our first check should compare the differences between the two computation methods. Both the branch and a local copy of develop were modified to dump the lat/lons to a netcdf file.

GeorgeGayno-NOAA commented 2 years ago

Check lat/lon differences between develop at 31271f7 and the branch at 9f66174.

Test 1 - Check the lat/lon differences for 13km RAP data. This data is on a rotated lat/lon 'B' grid (grid definition template number - 32769).

   grid_template=32769:winds(grid):
        I am not an Arakawa E-grid.
        I am rotated but have no rotation angle.
        I am staggered. What am I?
        (953 x 834) units 1e-06 input WE:SN output WE:SN res 56
        lat0 -10.590603 lat-center 54.000000 dlat 121.813000
        lon0 220.914154 lon-center 254.000000 dlon 121.813000 #points=794802

The wgrib2 library can't compute the lat/lons for this grid, so they are read into chgres_cube from a file - latlon_grid3.32769.nc.

The differences between the file and gdswzd computations are:

Variable        Group  Count         Sum      AbsSum          Min         Max      Range         Mean      StdDev
gridlat         /     790759    -120.266     162.903 -0.000530243 0.000289917 0.00082016 -0.000152089 0.000204953
gridlon         /     792752       49.76     771.789     -79.9549         360    439.955  6.27686e-05    0.414187
gridlat_corners /     792644    -120.655     164.868 -0.000541687  0.00219727 0.00273895 -0.000152219 0.000206831
gridlon_corners /     794497    -229.414     323.772    -0.181824    0.178787   0.360611 -0.000288754 0.000780303

The mean are very small compared to the grid resolution of 13 km. Some large differences in longitude are near the dateline, where the sign of longitude switches.

GeorgeGayno-NOAA commented 2 years ago

Test 2 - Check the lat/lon differences for the 3km HRRR data. This data is on a lambert conformal grid:

  grid_template=30:winds(grid):
        Lambert Conformal: (1799 x 1059) input WE:SN output WE:SN res 8
        Lat1 21.138123 Lon1 237.280472 LoV 262.500000
        LatD 38.500000 Latin1 38.500000 Latin2 38.500000
        LatSP 0.000000 LonSP 0.000000
        North Pole (1799 x 1059) Dx 3000.000000 m Dy 3000.000000 m mode 8

The differences between the wgrib2 and gdswzd computations are:

Variable                  Group   Count      Sum  AbsSum          Min          Max       Range         Mean      StdDev
gridlat                   /     1885229   98.535 149.762 -0.000123978  0.000217438 0.000341415  5.22669e-05  7.8939e-05
gridlon                   /     1905141 -871.205 871.205 -0.000732422 -0.000244141 0.000488281 -0.000457291 0.000110639
gridlat_corners           /     1907294  830.204 4353.38  -0.00643921    0.0065918    0.013031  0.000435279   0.0026419
gridlon_corners           /     1908000 -65904.6 65961.1   -0.0474243    0.0458069   0.0932312   -0.0345412  0.00527456

The mean differences at the center of the grid cell are very small compared to the grid resolution of 3 km. And the corner point latitude differences are small. However, the mean corner point longitude differences are comparable to the grid resolution. Something is not correct.

GeorgeGayno-NOAA commented 2 years ago

Test 3 - Check the lat/lon differences for 12km NAM data. This data is on a lambert conformal grid:

    grid_template=30:winds(grid):
        Lambert Conformal: (614 x 428) input WE:SN output WE:SN res 56
        Lat1 12.190000 Lon1 226.541000 LoV 265.000000
        LatD 25.000000 Latin1 25.000000 Latin2 25.000000
        LatSP -90.000000 LonSP 0.000000
        North Pole (614 x 428) Dx 12191.000000 m Dy 12191.000000 m mode 56

The differences between the wgrib2 and gdswzd computations are:

Variable                  Group  Count      Sum  AbsSum          Min         Max       Range        Mean      StdDev
gridlat                   /     261655  22.6143 23.3464 -6.58035e-05 0.000183105 0.000248909 8.64278e-05 5.58303e-05
gridlon                   /     254008  48.8392 50.8008 -0.000137329 0.000549316 0.000686646 0.000192274 0.000141188
gridlat_corners           /     263804  513.784 2193.29   -0.0313454   0.0270767   0.0584221   0.0019476  0.00970519
gridlon_corners           /     263835 -36764.1 36866.6    -0.212463    0.209259    0.421722   -0.139345   0.0287544

The mean differences at the center of the grid cell are very small compared to the grid resolution of 12 km. And the corner point latitude differences are small. However, as with the 3 km HRRR data, the mean corner point longitude differences are comparable to the grid resolution. Something is not correct.

GeorgeGayno-NOAA commented 2 years ago

Test 2 - Check the lat/lon differences for the 3km HRRR data. This data is on a lambert conformal grid:

  grid_template=30:winds(grid):
        Lambert Conformal: (1799 x 1059) input WE:SN output WE:SN res 8
        Lat1 21.138123 Lon1 237.280472 LoV 262.500000
        LatD 38.500000 Latin1 38.500000 Latin2 38.500000
        LatSP 0.000000 LonSP 0.000000
        North Pole (1799 x 1059) Dx 3000.000000 m Dy 3000.000000 m mode 8

The differences between the wgrib2 and gdswzd computations are:

Variable                  Group   Count      Sum  AbsSum          Min          Max       Range         Mean      StdDev
gridlat                   /     1885229   98.535 149.762 -0.000123978  0.000217438 0.000341415  5.22669e-05  7.8939e-05
gridlon                   /     1905141 -871.205 871.205 -0.000732422 -0.000244141 0.000488281 -0.000457291 0.000110639
gridlat_corners           /     1907294  830.204 4353.38  -0.00643921    0.0065918    0.013031  0.000435279   0.0026419
gridlon_corners           /     1908000 -65904.6 65961.1   -0.0474243    0.0458069   0.0932312   -0.0345412  0.00527456

The mean differences at the center of the grid cell are very small compared to the grid resolution of 3 km. And the corner point latitude differences are small. However, the mean corner point longitude differences are comparable to the grid resolution. Something is not correct.

Attaching some plots of corner point longitude minus center point longitude. The wgrib2 plot shows this difference is of the wrong sign (It is positive everywhere. For this N Amer. grid, it should be negative). The gdswzd plots show the expected magnitude (about 0.015 deg or 1.5 km) and sign.

hrrr wgrib2

hrrr gdswzd

GeorgeGayno-NOAA commented 2 years ago

Test 3 - Check the lat/lon differences for 12km NAM data. This data is on a lambert conformal grid:

    grid_template=30:winds(grid):
        Lambert Conformal: (614 x 428) input WE:SN output WE:SN res 56
        Lat1 12.190000 Lon1 226.541000 LoV 265.000000
        LatD 25.000000 Latin1 25.000000 Latin2 25.000000
        LatSP -90.000000 LonSP 0.000000
        North Pole (614 x 428) Dx 12191.000000 m Dy 12191.000000 m mode 56

The differences between the wgrib2 and gdswzd computations are:

Variable                  Group  Count      Sum  AbsSum          Min         Max       Range        Mean      StdDev
gridlat                   /     261655  22.6143 23.3464 -6.58035e-05 0.000183105 0.000248909 8.64278e-05 5.58303e-05
gridlon                   /     254008  48.8392 50.8008 -0.000137329 0.000549316 0.000686646 0.000192274 0.000141188
gridlat_corners           /     263804  513.784 2193.29   -0.0313454   0.0270767   0.0584221   0.0019476  0.00970519
gridlon_corners           /     263835 -36764.1 36866.6    -0.212463    0.209259    0.421722   -0.139345   0.0287544

The mean differences at the center of the grid cell are very small compared to the grid resolution of 12 km. And the corner point latitude differences are small. However, as with the 3 km HRRR data, the mean corner point longitude differences are comparable to the grid resolution. Something is not correct.

Some plots. As in the previous comment (https://github.com/ufs-community/UFS_UTILS/issues/591#issuecomment-1083477629), note how the gdswzd computation make more sense:

nam wgrib2

nam gdswzd

GeorgeGayno-NOAA commented 2 years ago

All files containing the input gridpoint latitude/longitudes and their differences (from the above tests) are on Hera here: /scratch1/NCEPDEV/da/George.Gayno/ufs_utils.git/chgres_cube_nowgrib2/PR641/latlon_check

Under each sub-directory, there are these netcdf files:

GeorgeGayno-NOAA commented 2 years ago

Graphics from all failed consistency tests are on Hera here: /scratch1/NCEPDEV/da/George.Gayno/ufs_utils.git/chgres_cube_nowgrib2/PR641

Under each test sub-directory, the graphics are here:

For this latest PR (#641) all differences are due solely to the difference between how the branch and develop compute the input grid latitude/longitudes.