NCAR / pynio

PyNIO is a multi-format data I/O package with a NetCDF-style interface
http://www.pyngl.ucar.edu/Nio.shtml
Apache License 2.0
112 stars 37 forks source link

forecast time not always read correctly from grib2 files #27

Open chiaral opened 5 years ago

chiaral commented 5 years ago

The grib files belong to SubX models, in particular, I refer to one of these models:

EMC/GEFS hindcast - pr_sfc_GEFS_06oct1999_00z_d01_d35_m00.grb

if I explore the data with wgrib2

wgrib2 pr_sfc_GEFS_06oct1999_00z_d01_d35_m00.grb2
1:0:d=1999100600:APCP:surface:1 day fcst:ENS=low-res ctl
2:71544:d=1999100600:APCP:surface:2 day fcst:ENS=low-res ctl
3:142213:d=1999100600:APCP:surface:3 day fcst:ENS=low-res ctl
4:211380:d=1999100600:APCP:surface:4 day fcst:ENS=low-res ctl
5:280114:d=1999100600:APCP:surface:5 day fcst:ENS=low-res ctl
6:348933:d=1999100600:APCP:surface:6 day fcst:ENS=low-res ctl

and so on. Note, it is a daily forecast. using wgrib2 with option -V


wgrib2 pr_sfc_GEFS_06oct1999_00z_d01_d35_m00.grb2 -V
1:0:vt=**1999100700**:surface:1 day fcst:APCP Total Precipitation [kg/m^2]:ENS=low-res ctl
    ndata=65160:undef=0:mean=2.85094:min=0:max=305.2
    grid_template=0:winds(N/S):
    lat-lon grid:(360 x 181) units 1e-06 input WE:NS output WE:SN res 48
    lat 90.000000 to -90.000000 by 1.000000
    lon 0.000000 to 359.000000 by 1.000000 #points=65160

2:71544:vt=**1999100800**:surface:2 day fcst:APCP Total Precipitation [kg/m^2]:ENS=low-res ctl
    ndata=65160:undef=0:mean=2.57202:min=0:max=320
    grid_template=0:winds(N/S):
    lat-lon grid:(360 x 181) units 1e-06 input WE:NS output WE:SN res 48
    lat 90.000000 to -90.000000 by 1.000000
    lon 0.000000 to 359.000000 by 1.000000 #points=65160

3:142213:vt=**1999100900**:surface:3 day fcst:APCP Total Precipitation [kg/m^2]:ENS=low-res ctl
    ndata=65160:undef=0:mean=2.51068:min=0:max=295.1
    grid_template=0:winds(N/S):
    lat-lon grid:(360 x 181) units 1e-06 input WE:NS output WE:SN res 48
    lat 90.000000 to -90.000000 by 1.000000
    lon 0.000000 to 359.000000 by 1.000000 #points=65160

and so on. Note the "vt" key, which is the validation time (? maybe) has the full calendar date of the lead time date.

When I open the same file both using Nio or xarray and the pynio engine I get:

Nio file:   pr_sfc_GEFS_06oct1999_00z_d01_d35_m00.grb2
   global attributes:
   dimensions:
      forecast_time0 = 35
      lat_0 = 181
      lon_0 = 360
   variables:
      float APCP_P1_L1_GLL0 [ forecast_time0, lat_0, lon_0 ]
         center :   US National Weather Service - NCEP (WMC)
         production_status :    Operational products
         long_name :    Total precipitation
         units :    kg m-2
         _FillValue :   1e+20
         grid_type :    Latitude/longitude
         parameter_discipline_and_category :    Meteorological products, Moisture
         parameter_template_discipline_category_number :    [1, 0, 1, 8]
         level_type :   Ground or water surface
         level :    0
         initial_time : 10/06/1999 (00:00)
      float lat_0 [ lat_0 ]
         long_name :    latitude
         grid_type :    Latitude/Longitude
         units :    degrees_north
         Dj :   1
         Di :   1
         Lo2 :  359
         La2 :  -90
         Lo1 :  0
         La1 :  90
      float lon_0 [ lon_0 ]
         long_name :    longitude
         grid_type :    Latitude/Longitude
         units :    degrees_east
         Dj :   1
         Di :   1
         Lo2 :  359
         La2 :  -90
         Lo1 :  0
         La1 :  90
      **integer forecast_time0 [ forecast_time0 ]
         long_name :    Forecast offset from initial time
         units :    hours**

and

<xarray.Dataset>
Dimensions:          (forecast_time0: 35, lat_0: 181, lon_0: 360)
Coordinates:
  * lat_0            (lat_0) float32 90.0 89.0 88.0 87.0 ... -88.0 -89.0 -90.0
  * lon_0            (lon_0) float32 0.0 1.0 2.0 3.0 ... 356.0 357.0 358.0 359.0
  * forecast_time0   (forecast_time0) timedelta64[ns] 01:00:00 ... 1 days 11:00:00
Data variables:
    APCP_P1_L1_GLL0  (forecast_time0, lat_0, lon_0) float32 ...

more specifically, it reads the forecast time as hours

temp.forecast_time0/np.timedelta64(1,'h')
<xarray.DataArray 'forecast_time0' (forecast_time0: 35)>
array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12., 13., 14.,
       15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.,
       29., 30., 31., 32., 33., 34., 35.])
Coordinates:
  * forecast_time0  (forecast_time0) timedelta64[ns] 01:00:00 ... 1 days 11:00:00

the same happens for the EMC/GEFS forecast data - so the same issue is in both the hindcast and realtime forecast data.

Any hint on how to fix this?