1

I have a GRIB dataset representing temperature over North America, in Lambert Conformal Conic (with specified central meridian, latitude of origin, etc.). It has what should be understood as a temporal dimension, given in Unix time. I am trying to use GDAL to both reproject this to EPSG:4326 and convert the file to NetCDF. However when I attempt to convert to NetCDF, the temporal dimension is not understood; instead each time dimension is interpreted as another band. This is leading to problems reading the output NetCDF in existing software, which expects to be able to look at the dimensions of a NetCDF header and then infer how the data is structured. How can I correctly use GDAL to perform this conversion, with the output being an appropriately-structured NetCDF with temporal dimension?

Pruned output of gdalinfo input_dataset.bin:

Driver: GRIB/GRIdded Binary (.grb)
Files: ds.temp.bin
Size is 2145, 1377
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6371200,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",25],
    PARAMETER["standard_parallel_2",25],
    PARAMETER["latitude_of_origin",25],
    PARAMETER["central_meridian",265],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
PROJ.4 string is:
'+proj=lcc +lat_1=25 +lat_2=25 +lat_0=25 +lon_0=265 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs '
Origin = (-2764474.350731994491071,3232111.710792394354939)
Pixel Size = (2539.702999999999975,-2539.702999999999975)
Corner Coordinates:
Upper Left  (-2764474.351, 3232111.711) (130d 7'22.60"W, 49d56'49.26"N)
Lower Left  (-2764474.351, -265059.320) (121d33'48.72"W, 20d10'43.04"N)
Upper Right ( 2683188.584, 3232111.711) ( 60d51'57.95"W, 50d 6'46.52"N)
Lower Right ( 2683188.584, -265059.320) ( 69d11'54.73"W, 20d19' 6.40"N)
Center      (  -40642.883, 1483526.195) ( 95d27' 8.65"W, 38d13' 5.87"N)
Band 1 Block=2145x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC="Ground or water surface"
  NoData Value=9999
  Metadata:
    GRIB_COMMENT=Temperature [C]
    GRIB_ELEMENT=T
    GRIB_FORECAST_SECONDS=3600 sec
    GRIB_REF_TIME=  1455177600 sec UTC
    GRIB_SHORT_NAME=0-SFC
    GRIB_UNIT=[C]
    GRIB_VALID_TIME=  1455181200 sec UTC
Band 2 Block=2145x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC="Ground or water surface"
  NoData Value=9999
  Metadata:
    GRIB_COMMENT=Temperature [C]
    GRIB_ELEMENT=T
    GRIB_FORECAST_SECONDS=7200 sec
    GRIB_REF_TIME=  1455177600 sec UTC
    GRIB_SHORT_NAME=0-SFC
    GRIB_UNIT=[C]
    GRIB_VALID_TIME=  1455184800 sec UTC
Band 3 Block=2145x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC="Ground or water surface"
  NoData Value=9999
  Metadata:
    GRIB_COMMENT=Temperature [C]
    GRIB_ELEMENT=T
    GRIB_FORECAST_SECONDS=10800 sec
    GRIB_REF_TIME=  1455177600 sec UTC
    GRIB_SHORT_NAME=0-SFC
    GRIB_UNIT=[C]
    GRIB_VALID_TIME=  1455188400 sec UTC

gdal_warp command:

gdalwarp -s_srs '+proj=lcc +lat_1=25 +lat_2=25 +lat_0=25 +lon_0=265 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs ' -t_srs EPSG:4326 -of netCDF input_dataset.bin output_dataset.nc

Output of gdalinfo output_dataset.nc

Driver: netCDF/Network Common Data Format
Files: output_dataset.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#Conventions=CF-1.5
  NC_GLOBAL#GDAL=GDAL 1.11.2, released 2015/02/10
  NC_GLOBAL#history=Mon Feb 15 17:07:52 2016: GDAL CreateCopy( temp.nc, ... )
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"temp.nc":Band1
  SUBDATASET_1_DESC=[1377x2145] Band1 (64-bit floating-point)
  SUBDATASET_2_NAME=NETCDF:"temp.nc":Band2
  SUBDATASET_2_DESC=[1377x2145] Band2 (64-bit floating-point)
  SUBDATASET_3_NAME=NETCDF:"temp.nc":Band3
  SUBDATASET_3_DESC=[1377x2145] Band3 (64-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

(Note, I'm also confused as to why the output has no specified coordinate system.)

In the past I have used wgrib for this conversion, but I have encountered problems with reprojection so I'd like to work out the equivalent GDAL command. A desired output structure (this from a semi-successful wgrib conversion) would be:

Driver: netCDF/Network Common Data Format
Files: outupt_dataset.nc
Size is 2651, 1551
Coordinate System is `'
Origin = (-122.010000000000005,51.009999999999998)
Pixel Size = (0.020000000000010,-0.020000000000000)
Metadata:
  latitude#long_name=latitude
  latitude#units=degrees_north
  longitude#long_name=longitude
  longitude#units=degrees_east
  NC_GLOBAL#Conventions=COARDS
  NC_GLOBAL#GRIB2_grid_template=0
  NC_GLOBAL#History=created by wgrib2
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={46,6}
  NETCDF_DIM_time_VALUES={1455181200,1455184800,1455188400,1455192000,1455195600,1455199200,1455202800,1455206400,1455210000,1455213600,1455217200,1455220800,1455224400,1455228000,1455231600,1455235200,1455238800,1455242400,1455246000,1455249600,1455253200,1455256800,1455260400,1455264000,1455267600,1455271200,1455274800,1455278400,1455282000,1455285600,1455289200,1455292800,1455296400,1455300000,1455303600,1455307200,1455310800,1455321600,1455332400,1455343200,1455354000,1455364800,1455375600,1455386400,1455397200,1455408000}
  time#long_name=verification time generated by wgrib2 function verftime()
  time#reference_date=2016.02.11 08:00:00 UTC
  time#reference_time=1455177600
  time#reference_time_description=forecast or accumulated, reference date is fixed
  time#reference_time_type=3
  time#time_step=-1
  time#time_step_setting=auto
  time#units=seconds since 1970-01-01 00:00:00.0 0:00
  TMP_surface#_FillValue=9.9990003e+20
  TMP_surface#level=surface
  TMP_surface#long_name=Temperature
  TMP_surface#short_name=TMP_surface
  TMP_surface#units=K
Corner Coordinates:
Upper Left  (-122.0100000,  51.0100000) 
Lower Left  (-122.0100000,  19.9900000) 
Upper Right ( -68.9900000,  51.0100000) 
Lower Right ( -68.9900000,  19.9900000) 
Center      ( -95.5000000,  35.5000000) 
Band 1 Block=2651x1551 Type=Float32, ColorInterp=Undefined
  NoData Value=9.99900026055400882e+20
  Metadata:
    _FillValue=9.9990003e+20
    level=surface
    long_name=Temperature
    NETCDF_DIM_time=1455181200
    NETCDF_VARNAME=TMP_surface
    short_name=TMP_surface
    units=K
Band 2 Block=2651x1551 Type=Float32, ColorInterp=Undefined
  NoData Value=9.99900026055400882e+20
  Metadata:
    _FillValue=9.9990003e+20
    level=surface
    long_name=Temperature
    NETCDF_DIM_time=1455184800
    NETCDF_VARNAME=TMP_surface
    short_name=TMP_surface
    units=K
Band 3 Block=2651x1551 Type=Float32, ColorInterp=Undefined
  NoData Value=9.99900026055400882e+20
  Metadata:
    _FillValue=9.9990003e+20
    level=surface
    long_name=Temperature
    NETCDF_DIM_time=1455188400
    NETCDF_VARNAME=TMP_surface
    short_name=TMP_surface
    units=K

Notice in the last header, there is this dimension information:

NETCDF_DIM_EXTRA={time}
NETCDF_DIM_time_DEF={46,6}

I want gdalwarp or gdal_translate to produce something equivalent in the dimensions (not using bands). Do I set these after the fact with gdal_translate, or is there a more automated way to determine the temporal dimensions?

A related question has been asked previously (Converting GRIB to NetCDF with time dimension using netcdfAll?) but that is about using the netcdfAll Java library, when I want to use GDAL, if possible.

Improve this question
2
  • I don't believe you can create 3D variables with GDAL - see "creation issues": gdal.org/frmt_netcdf.html You'll have to recompose them with something like the nco operators or R/python etc. or convert them from Grib some other way. (I would use R, but then creating NetCDF files makes me feel terrible so I try to avoid it).
    – mdsumner
    Commented Feb 15, 2016 at 8:50
  • Huh. Until this very moment I thought GDAL was all powerful. I'll accept that as an answer if you write something up, given what I just read in that documentation ("Each variable array is named Band1, Band2, ... BandN.").
    – alphabetasoup
    Commented Feb 15, 2016 at 8:59

1 Answer 1

Reset to default
1

I don't believe you can create 3D variables with GDAL - see "creation issues" here:

http://www.gdal.org/frmt_netcdf.html

You'll have to recompose them with something like the nco operators or R/python etc. or convert them from Grib some other way.

FWIW, it's about a 2 minute job in R as long as you can read and make sense of the grib file: 

library(raster)
library(rgdal)  ## will need Grib available see gdalDrivers()

## hope for the best
b <- brick("input_dataset.bin")

## projectRaster(b, crs = "+init=EPSG:4326"), other trickery here 
## ... etc.
##

writeRaster(b, "output_dataset.nc")

Something like that will give you a useable solution, but it really depends on the details for a particular data set. Can you point to an example file?

Please note that converting from a projection to longlat is probably needlessly destructive, why are you smashing this nice map projection?

Improve this answer
1
  • No worries, I have other ways to convert, I just find GDAL typically the easiest option. Regarding the projection, at some point the data will be be projected to EPSG:900913: the NetCDF is input data for a tile server. My other datasets are EPSG:4327, I just wanted to avoid improving my projection logic.
    – alphabetasoup
    Commented Feb 15, 2016 at 9:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.