2

I've prepared a netcdf4 file called dummy_dataset.nc, using the netcdf4 python library:

from netCDF4    import Dataset, num2date, date2num
from numpy      import linspace, zeros, random
from datetime import datetime, timedelta

ds   = Dataset("dummy_dataset_2.nc", "w", clobber = True)
zoom = 4
coef = 2**zoom
xsize=360*coef
ysize=180*coef
tsize=1
lsize=1
ds.createDimension("lat", ysize)  
ds.createDimension("lon", xsize)
ds.createDimension("time", tsize)
ds.createDimension("level", lsize)

latitudes   = ds.createVariable("lat", "f4", ('y',))
longitudes  = ds.createVariable("lon", "f4", ('x',))
times       = ds.createVariable("time","f8", ('time',))
levels      = ds.createVariable("level", "i4", ('level',))


latitudes.units  = "degrees_north"
longitudes.units = "degrees_east"
longitudes.axis = "x"
latitudes.axis = "y"
times.units      = "days since 1860-01-01 00:00:00"
times.calendar   = "standard"
levels.units = "Km"

dates = [datetime(2001,3,1)+n*timedelta(hours=12) for n in range(tsize)]

latitudes[:]    = linspace(-90,90,ysize, endpoint=True)
longitudes[:]   = linspace(-180,180, xsize, endpoint=True)
times[:]        = date2num(dates,times.units,calendar=times.calendar)
levels[:]       = range(lsize)


var1 = ds.createVariable("test0", "f4", ("lat","lon",),chunksizes=(coef,coef,))
var3 = ds.createVariable("test1", "f4", ('time','level','lat','lon',),chunksizes=(1,1,coef,coef,))

var1[:] = random.random((ysize,xsize))
var3[:] = random.random((tsize,lsize,ysize,xsize))

ds.Conventions = "CF-1.0"

ds.close()

This file has to be processed with gdal functions, written using the python bindings. The problem is that it doesn't seem to get the geotransform correctly. In fact when i try to "gdalinfo" a file subset(let's say test0): 

Driver: netCDF/Network Common Data Format
Files: dummy_dataset.nc
       dummy_dataset.nc.aux.xml
Size is 5760, 2880
Coordinate System is `'
Metadata:
  level#units=Km
  NC_GLOBAL#Conventions=CF-1.0
  time#calendar=standard
  time#units=days since 1860-01-01 00:00:00
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 2880.0)
Upper Right ( 5760.0,    0.0)
Lower Right ( 5760.0, 2880.0)
Center      ( 2880.0, 1440.0)
Band 1 Block=5760x1 Type=Float32, ColorInterp=Undefined
  Min=0.000 Max=1.000 
  Minimum=0.000, Maximum=1.000, Mean=0.500, StdDev=0.289
  NoData Value=9.96920996838686905e+36
  Metadata:
    NETCDF_VARNAME=test0
    STATISTICS_MAXIMUM=0.99999994039536
    STATISTICS_MEAN=0.50011895194454
    STATISTICS_MINIMUM=2.25937498044e-07
    STATISTICS_STDDEV=0.28865938167717



Files: dummy_dataset.nc
Size is 5760, 2880
Coordinate System is `'
Metadata:
  level#units=Km
  NC_GLOBAL#Conventions=CF-1.0
  time#calendar=standard
  time#units=days since 1860-01-01 00:00:00
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 2880.0)
Upper Right ( 5760.0,    0.0)
Lower Right ( 5760.0, 2880.0)
Center      ( 2880.0, 1440.0)
Band 1 Block=5760x1 Type=Float32, ColorInterp=Undefined
  Minimum=0.000, Maximum=1.000, Mean=0.500, StdDev=0.289
  NoData Value=9.96920996838686905e+36
  Metadata:
    NETCDF_VARNAME=test0
    STATISTICS_MAXIMUM=0.99999994039536
    STATISTICS_MEAN=0.50011895194454
    STATISTICS_MINIMUM=2.25937498044e-07
    STATISTICS_STDDEV=0.28865938167717

noticeably, the size of the grid is printed instead of the desired coverage in terms of bounding box (in this case lon[-180, 180],lat[-90, 90]). If we print the geotransform of the file, we get the following:

---NETCDF:"dummy_dataset.nc":test0---
Driver: netCDF/Network Common Data Format
Size is 5760, 2880
GetGeoTransform() =  (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
---NETCDF:"dummy_dataset.nc":test1---
Driver: netCDF/Network Common Data Format
Size is 5760, 2880
GetGeoTransform() =  (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)

At first, i thought it was a driver limitation on gdal part. So i downloaded the classic ECMWF_ERA-40_subset.nc and tested against it.

Here's the result:

Driver: netCDF/Network Common Data Format
Files: ../files/ECMWF_ERA-40_subset.nc
Size is 512, 512
Coordinate System is `'
Metadata:
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#history=2004-09-15 17:04:29 GMT by mars2netcdf-0.92
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":tcw
  SUBDATASET_1_DESC=[62x73x144] tcw (16-bit integer)
  SUBDATASET_2_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":tcwv
  SUBDATASET_2_DESC=[62x73x144] tcwv (16-bit integer)
  SUBDATASET_3_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":lsp
  SUBDATASET_3_DESC=[62x73x144] lsp (16-bit integer)
  SUBDATASET_4_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":cp
  SUBDATASET_4_DESC=[62x73x144] cp (16-bit integer)
  SUBDATASET_5_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":msl
  SUBDATASET_5_DESC=[62x73x144] msl (16-bit integer)
  SUBDATASET_6_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":blh
  SUBDATASET_6_DESC=[62x73x144] blh (16-bit integer)
  SUBDATASET_7_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":tcc
  SUBDATASET_7_DESC=[62x73x144] tcc (16-bit integer)
  SUBDATASET_8_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":p10u
  SUBDATASET_8_DESC=[62x73x144] p10u (16-bit integer)
  SUBDATASET_9_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":p10v
  SUBDATASET_9_DESC=[62x73x144] p10v (16-bit integer)
  SUBDATASET_10_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":p2t
  SUBDATASET_10_DESC=[62x73x144] p2t (16-bit integer)
  SUBDATASET_11_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":p2d
  SUBDATASET_11_DESC=[62x73x144] p2d (16-bit integer)
  SUBDATASET_12_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":e
  SUBDATASET_12_DESC=[62x73x144] e (16-bit integer)
  SUBDATASET_13_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":lcc
  SUBDATASET_13_DESC=[62x73x144] lcc (16-bit integer)
  SUBDATASET_14_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":mcc
  SUBDATASET_14_DESC=[62x73x144] mcc (16-bit integer)
  SUBDATASET_15_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":hcc
  SUBDATASET_15_DESC=[62x73x144] hcc (16-bit integer)
  SUBDATASET_16_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":tco3
  SUBDATASET_16_DESC=[62x73x144] tco3 (16-bit integer)
  SUBDATASET_17_NAME=NETCDF:"../files/ECMWF_ERA-40_subset.nc":tp
  SUBDATASET_17_DESC=[62x73x144] tp (16-bit integer)
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)

and here for the subdataset:

Driver: netCDF/Network Common Data Format
Files: ../files/ECMWF_ERA-40_subset.nc
Size is 144, 73
Coordinate System is `'
Origin = (-1.250000000000000,91.250000000000000)
Pixel Size = (2.500000000000000,-2.500000000000000)
Metadata:
  latitude#long_name=latitude
  latitude#units=degrees_north
  longitude#long_name=longitude
  longitude#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.0
  NC_GLOBAL#history=2004-09-15 17:04:29 GMT by mars2netcdf-0.92
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={62,4}
  NETCDF_DIM_time_VALUES={898476,898482,898500,898506,898524,898530,898548,898554,898572,898578,898596,898602,898620,898626,898644,898650,898668,898674,898692,898698,898716,898722,898740,898746,898764,898770,898788,898794,898812,898818,898836,898842,898860,898866,898884,898890,898908,898914,898932,898938,898956,898962,898980,898986,899004,899010,899028,899034,899052,899058,899076,899082,899100,899106,899124,899130,899148,899154,899172,899178,899196,899202}
  tcw#add_offset=44.3250482744756
  tcw#long_name=Total column water
  tcw#missing_value=-32767
  tcw#scale_factor=0.001350098174548095
  tcw#units=kg m**-2
  tcw#_FillValue=-32767
  time#long_name=time
  time#units=hours since 1900-01-01 00:00:0.0
Corner Coordinates:
Upper Left  (      -1.250,      91.250) 
Lower Left  (      -1.250,     -91.250) 
Upper Right (     358.750,      91.250) 
Lower Right (     358.750,     -91.250) 
Center      ( 178.7500000,   0.0000000) 
Band 1 Block=144x1 Type=Int16, ColorInterp=Undefined
  NoData Value=-32767
  Unit Type: kg m**-2
  Offset: 44.3250482744756,   Scale:0.0013500981745481
  Metadata:
    add_offset=44.3250482744756
    long_name=Total column water
    missing_value=-32767
    NETCDF_DIM_time=898476
    NETCDF_VARNAME=tcw
    scale_factor=0.001350098174548095
    units=kg m**-2
    _FillValue=-32767

Notice how in the case the subset extents are properly shown and actually computed.

---NETCDF:"../files/ECMWF_ERA-40_subset.nc":tcw---
Driver: netCDF/Network Common Data Format
Size is 144, 73
GetGeoTransform() =  (-1.25, 2.5, 0.0, 91.25, 0.0, -2.5)
---NETCDF:"../files/ECMWF_ERA-40_subset.nc":tcwv---
Driver: netCDF/Network Common Data Format
Size is 144, 73
GetGeoTransform() =  (-1.25, 2.5, 0.0, 91.25, 0.0, -2.5)
---NETCDF:"../files/ECMWF_ERA-40_subset.nc":lsp---
Driver: netCDF/Network Common Data Format
Size is 144, 73
GetGeoTransform() =  (-1.25, 2.5, 0.0, 91.25, 0.0, -2.5)
---NETCDF:"../files/ECMWF_ERA-40_subset.nc":cp---
Driver: netCDF/Network Common Data Format
Size is 144, 73
GetGeoTransform() =  (-1.25, 2.5, 0.0, 91.25, 0.0, -2.5)
---NETCDF:"../files/ECMWF_ERA-40_subset.nc":msl---
Driver: netCDF/Network Common Data Format
Size is 144, 73
GetGeoTransform() =  (-1.25, 2.5, 0.0, 91.25, 0.0, -2.5)
...

now there are two interesting differences between the file i've produced and the one in analysis. The netcdfdump command proves to be helpful in this regard:

netcdf ECMWF_ERA-40_subset {
dimensions:
    longitude = 144 ;
    latitude = 73 ;
    time = UNLIMITED ; // (62 currently)
variables:
    float longitude(longitude) ;
        longitude:units = "degrees_east" ;
        longitude:long_name = "longitude" ;
    float latitude(latitude) ;
        latitude:units = "degrees_north" ;
        latitude:long_name = "latitude" ;
    int time(time) ;
        time:units = "hours since 1900-01-01 00:00:0.0" ;
        time:long_name = "time" ;
    short tcw(time, latitude, longitude) ;
        tcw:scale_factor = 0.0013500981745481 ;
        tcw:add_offset = 44.3250482744756 ;
        tcw:_FillValue = -32767s ;
        tcw:missing_value = -32767s ;
        tcw:units = "kg m**-2" ;
        tcw:long_name = "Total column water" ;
    short tcwv(time, latitude, longitude) ;
        tcwv:scale_factor = 0.001327110772669 ;
        tcwv:add_offset = 43.5704635546154 ;
        tcwv:_FillValue = -32767s ;
        tcwv:missing_value = -32767s ;
        tcwv:units = "kg m**-2" ;
        tcwv:long_name = "Total column water vapour" ;
    short lsp(time, latitude, longitude) ;
        lsp:scale_factor = 8.03329303850659e-07 ;
        lsp:add_offset = 0.0263210846406669 ;
        lsp:_FillValue = -32767s ;
        lsp:missing_value = -32767s ;
        lsp:units = "m" ;
        lsp:long_name = "Stratiform precipitation (Large-scale precipitation)" ;
  • My file is a netcdf4 and ECMWF_ERA is a netcdf classic
  • the presence of two metadata attributes: scale_factor and add_offset

The second one is probably the most relevant, however i have no clue on how to determine those factors to describe the kind of grid used in the file. Even once computed, are there some other "hidden" factors that affect the description (thus the understanding on gdal behalf) of the grid/grids saved in the file, and how should i determine their values? how can i represent (metadata wise), for example, a regular grid (such as in this case) with a fixed increment in degrees along both axis ?

Improve this question
2
  • Running the first bit of your code I get the following error ValueError("cannot find dimension %s in this group or parent groups" % dimname) You should change these two lines to latitudes = ds.createVariable("lat", "f4", ('lat',)) longitudes = ds.createVariable("lon", "f4", ('lon',)) Loading the test0 with gdal: import gdal ds = gdal.Open('NETCDF:dummy_dataset_2.nc:test0') the geot = ds_ph.GetGeoTransform() returns >>(-180.03125542628928, 0.06251085257857267, 0.0, 90.03126085446335, 0.0, -0.06252170892671066) which seems correct.
    – Andrea Massetti
    Commented Mar 7, 2018 at 0:06
  • @AndreaMassetti you were absolutely right. I forgot that dimensions must have the same name as their associated variables. However something is still boggling me.
    – darkpirate
    Commented Mar 8, 2018 at 13:37

1 Answer 1

Reset to default
2

As AndreaMassetti pointed out, the code i presented was malformed, i was messing around with variables naming to check if it would actually make a difference, and in a sense, it does. Not much in the name of the variable itself but in the fact that to reflect the C.F. naming convention, dimension and their associated variable must have the same name.

netcdf dummy_dataset_2 {
dimensions:
    lat = 360 ;
    lon = 720 ;
    time = 2 ;
    level = 2 ;
variables:
    double lat(lat) ;
        lat:units = "degrees_north" ;
        lat:axis = "y" ;
        lat:standard_name = "latitude" ;
    double lon(lon) ;
        lon:units = "degrees_east" ;
        lon:axis = "x" ;
        lon:standard_name = "longitude" ;
    double time(time) ;
        time:units = "days since 1860-01-01 00:00:00" ;
        time:calendar = "standard" ;
    int level(level) ;
        level:units = "Km" ;
    double test_temp(time, level, lat, lon) ;
        test_temp:units = "K" ;
        test_temp:standard_name = "temperature" ;

// global attributes:
        :Conventions = "CF-1.6" ;
}

and the resulting gdal info:

Files: dummy_dataset_2.nc
Size is 720, 360
Coordinate System is `'
Origin = (-180.250000000000000,90.250000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  lat#axis=y
  lat#standard_name=latitude
  lat#units=degrees_north
  level#units=Km
  lon#axis=x
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.6
  NETCDF_DIM_EXTRA={time,level}
  NETCDF_DIM_level_DEF={2,4}
  NETCDF_DIM_level_VALUES={0,1}
  NETCDF_DIM_time_DEF={2,6}
  NETCDF_DIM_time_VALUES={51559,51559.5}
  test_temp#standard_name=temperature
  test_temp#units=K
  time#calendar=standard
  time#units=days since 1860-01-01 00:00:00
Corner Coordinates:
Upper Left  (-180.2500000,  90.2500000) 
Lower Left  (-180.2500000, -89.7500000) 
Upper Right ( 179.7500000,  90.2500000) 
Lower Right ( 179.7500000, -89.7500000) 
Center      (  -0.2500000,   0.2500000) 
Band 1 Block=720x1 Type=Float64, ColorInterp=Undefined
  NoData Value=9.96920996838686905e+36
  Unit Type: K
  Metadata:
    NETCDF_DIM_level=0
    NETCDF_DIM_time=51559
    NETCDF_VARNAME=test_temp
    standard_name=temperature
    units=K
Band 2 Block=720x1 Type=Float64, ColorInterp=Undefined
  NoData Value=9.96920996838686905e+36
  Unit Type: K
  Metadata:
    NETCDF_DIM_level=1
    NETCDF_DIM_time=51559
    NETCDF_VARNAME=test_temp
    standard_name=temperature
    units=K
Band 3 Block=720x1 Type=Float64, ColorInterp=Undefined
  NoData Value=9.96920996838686905e+36
  Unit Type: K
  Metadata:
    NETCDF_DIM_level=0
    NETCDF_DIM_time=51559.5
    NETCDF_VARNAME=test_temp
    standard_name=temperature
    units=K
Band 4 Block=720x1 Type=Float64, ColorInterp=Undefined
  NoData Value=9.96920996838686905e+36
  Unit Type: K
  Metadata:
    NETCDF_DIM_level=1
    NETCDF_DIM_time=51559.5
    NETCDF_VARNAME=test_temp
    standard_name=temperature
    units=K
Improve this answer

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.