ecmwf / cfgrib

A Python interface to map GRIB files to the NetCDF Common Data Model following the CF Convention using ecCodes
Apache License 2.0
400 stars 77 forks source link

GFS: Missing levels for sigmaLayer and depthBelowLandLayer #195

Closed maxhollmann closed 3 years ago

maxhollmann commented 3 years ago

cfgrib seems to miss some levels for the level types sigmaLayer and depthBelowLandLayer.

I'm working with this grib file: https://noaa-gfs-bdp-pds.s3.amazonaws.com/gfs.20210115/00/gfs.t00z.pgrb2.0p25.f003

wgrib2 -lev ~/Downloads/gfs.t00z.pgrb2.0p25.f003 shows four levels for each of these level types:

425:224793214:0-0.1 m below ground
427:225650288:0.1-0.4 m below ground
429:226508191:0.4-1 m below ground
431:227363072:1-2 m below ground

558:309354159:0.33-1 sigma layer
559:310143600:0.44-1 sigma layer
560:310922234:0.72-0.94 sigma layer
561:311736866:0.44-0.72 sigma layer

Reading the same file with cfgrib results in only two:

import cfgrib

ds = cfgrib.open_datasets("gfs.t00z.pgrb2.0p25.f003")

d = next(d for d in ds if 'depthBelowLandLayer' in d.variables)
d.depthBelowLandLayer
# <xarray.IndexVariable 'depthBelowLandLayer' (depthBelowLandLayer: 2)>
# array([0, 1])
# Attributes:
#     long_name:      soil depth
#     units:          m
#     positive:       down
#     standard_name:  depth

d = next(d for d in ds if 'sigmaLayer' in d.variables)
d.sigmaLayer
# <xarray.IndexVariable 'sigmaLayer' (sigmaLayer: 2)>
# array([0, 1])
# Attributes:
#     long_name:  original GRIB coordinate for key: level(sigmaLayer)
#     units:      1

d.r.shape
# (2, 721, 1440)
shahramn commented 3 years ago

Several of the messages in this file use non-standard i.e. local values for "typeOfFirstFixedSurface" and therefore the ecCodes key "typeOfLevel" returns "unknown".

% grib_ls -p count,typeOfFirstFixedSurface,typeOfSecondFixedSurface,typeOfLevel gfs.t00z.pgrb2.0p25.f003 |\
               grep unknown
8                         220                       255                       unknown
9                         220                       255                       unknown
10                        220                       255                       unknown
476                       200                       255                       unknown
477                       200                       255                       unknown
478                       200                       255                       unknown
479                       200                       255                       unknown
480                       214                       255                       unknown
548                       204                       255                       unknown

The key typeOfFirstFixedSurface and typeOfSecondFixedSurface use Code Table 4.5 in which we have:

...
184 184 Grid tile glacier ice and inland ice fraction as a model surface
# 185-191 Reserved 
# 192-254 Reserved for local use 
255 255 Missing 

So the values used in this file are in the range 192-254 i.e. they are local

maxhollmann commented 3 years ago

Right, but that's not what I mean here. The level types sigmaLayer and depthBelowLandLayer both get recognized. The values for these two levels are different when parsed by eccodes than by wgrib2:

grib_ls -p count,typeOfLevel,level gfs.t00z.pgrb2.0p25.f003 | grep -E "depthBelowLandLayer|sigmaLayer":

424          depthBelowLandLayer  0           
425          depthBelowLandLayer  0           
426          depthBelowLandLayer  0           
427          depthBelowLandLayer  0           
428          depthBelowLandLayer  0           
429          depthBelowLandLayer  0           
430          depthBelowLandLayer  1           
431          depthBelowLandLayer  1           
558          sigmaLayer   0           
559          sigmaLayer   0           
560          sigmaLayer   1           
561          sigmaLayer   0           

wgrib2 -lev gfs.t00z.pgrb2.0p25.f003 | grep -E "sigma layer|below ground":

424:224283344:0-0.1 m below ground
425:224793214:0-0.1 m below ground
426:225146662:0.1-0.4 m below ground
427:225650288:0.1-0.4 m below ground
428:226004347:0.4-1 m below ground
429:226508191:0.4-1 m below ground
430:226864070:1-2 m below ground
431:227363072:1-2 m below ground
558:309354159:0.33-1 sigma layer
559:310143600:0.44-1 sigma layer
560:310922234:0.72-0.94 sigma layer
561:311736866:0.44-0.72 sigma layer
shahramn commented 3 years ago

OK. Thanks for the clarification. You have to decode the "level" keys as DOUBLEs rather than INTEGERs Try:

% grib_ls -p count,typeOfLevel,topLevel:d,bottomLevel:d gfs.t00z.pgrb2.0p25.f003 |\
    grep -E "depthBelowLandLayer|sigmaLayer"
424          depthBelowLandLayer  0            0.1
425          depthBelowLandLayer  0            0.1
426          depthBelowLandLayer  0.1          0.4
427          depthBelowLandLayer  0.1          0.4
428          depthBelowLandLayer  0.4          1
429          depthBelowLandLayer  0.4          1
430          depthBelowLandLayer  1            2
431          depthBelowLandLayer  1            2
558          sigmaLayer   0.33         1
559          sigmaLayer   0.44         1
560          sigmaLayer   0.72         0.94
561          sigmaLayer   0.44         0.72

Here I am using the keys "topLevel" and "bottomLevel" decoded as doubles (using the ":d" qualifier)

alexamici commented 3 years ago

@shahramn then this looks like an actual bug (or missing feature) in cfgrib as it only requests the level key and it accepts the default data type, that is int.

Thanks for the clarification.

alexamici commented 3 years ago

Levels will be read as floats in the next release.