[netcdfgroup] how to output a netcdf with an irregular grid

• To: netcdfgroup@xxxxxxxxxxxxxxxx
• Subject: [netcdfgroup] how to output a netcdf with an irregular grid
• From: Michael Nolde <nolde@xxxxxxxxxxxxxxxxxxxxxx>
• Date: Thu, 27 Feb 2014 19:18:01 +0100

Dear list members,

I have been trying to write a netcdf with an irregular grid for weeks now,

the data in my output file to my input file (with correct irregular grid).
My current approch is to create a regular grid within Python and assign the

Here's the details:

I've got a netCDF file which I want to read, modify, and save, which works
well with R ('ncdf4') and Python ('python-netcdf').
The problem is, that the grid representing the data (temperature values) is
irregular, and I am not able to output an irregular grid. In the original
file I have two variables ('lat' and 'lon'), which are both two-dimensional
(356 x 236), with each entry representing the latitude or longitude,
respectively, of each cell in the original grid (356 x 236 x 365 [days]).

In the input file, the data grid (TMAX_2M, please see blow) references the
'lat' and 'lon' variables by the argument 'coordinates = "lon lat"', as far
as I understood, but I have been unable to produce the desired behaviour
even if I hand the datagrid with the 'coordinates'-flag to the output file.

I would love to know how to make the netCDF file aware that the
latitude/longitude coordinates of the data grid cells are stored in the
"external" 'lat' and 'lon' variables, and not in the grid object itself.

as well as my source code (in 'R') so far.
Any help is highly appreciated, thank you very much in advance.



Best regards, Michael


this is the structure of the source netCDF file:
-----------------------------------------------------

dimensions:
    x = 356;
    y = 236;
    height_2m = 1;
    time = UNLIMITED; // (365 currently
    nb2 = 2;

variables:
    float lon(y=236, x=356);
        :standard_name = "longitude";
        :long_name = "longitude";
        :units = "degrees_east";
        :_CoordinateAxisType = "Lon";

    float lat(y=236, x=356);
        :standard_name = "latitude";
        :long_name = "latitude";
        :units = "degrees_north";
        :_CoordinateAxisType = "Lat";

    float height_2m(height_2m=1);
        :standard_name = "height";
        :long_name = "height above the surface";
        :units = "m";
        :positive = "up";
        :axis = "Z";

    double time(time=365);
        :standard_name = "time";
        :long_name = "time";
        :bounds = "time_bnds";
        :units = "seconds since 1979-01-01 00:00:00";
        :calendar = "proleptic_gregorian";

    double time_bnds(time=365, nb2=2);
        :units = "seconds since 1979-01-01 00:00:00";
        :calendar = "proleptic_gregorian";

    float TMAX_2M(time=365, height_2m=1, y=236, x=356);
        :standard_name = "air_temperature";
        :long_name = "2m maximum temperature";
        :units = "K";
        :coordinates = "lon lat";
        :cell_methods = "time: maximum";
        :_FillValue = -9.0E33f; // float



this is my 'Python' - code so far,
based on an example by Jeff Whitaker (found on the Unidata-website)
----------------------------------------------------------------------------------------------
# the Scientific Python netCDF 3 interface
# http://dirac.cnrs-orleans.fr/ScientificPython/
from Scientific.IO.NetCDF import NetCDFFile as Dataset

# the 'classic' version of the netCDF4 python interface
# http://code.google.com/p/netcdf4-python/
#from netCDF4_classic import Dataset
from numpy import arange, dtype # array module from http://numpy.scipy.org


import netCDF4
import numpy as np
import sys

f = netCDF4.Dataset('TMAX_2M_daily1981_box.nc', 'r')
for v in f.variables:
    print(v),

type(f.variables) # f.variables is a dictionary data structure

for d in f.dimensions:
    print(d),
print


x = f.dimensions['x']
y = f.dimensions['y']
#time = f.dimensions['time']

temp = f.variables['TMAX_2M']
lon = f.variables['lon']
lat = f.variables['lat']
time = f.variables['time']
height = f.variables['height_2m']


npLon = np.array(lon)
npLat = np.array(lat)

npLonMin = npLon.min()
npLonMax = npLon.max()
npLatMin = npLat.min()
npLatMax = npLat.max()


lonOut = np.arange(-18.75, 59.00, 0.25)
latOut = np.arange(24.75, 60.25, 0.25)


tempOut = np.zeros(shape=(2,len(latOut),len(lonOut)),dtype='float32')
tempOut.fill(np.nan)



iTime = 0
iHeight = 0
tempOutOld = 0

for iY in range(0,236):
    for iX in range(0,356):
        cellX = int(round((lon[iY,iX] - npLonMin) * 4,0))
        cellY = int(round((lat[iY,iX] - npLatMin) * 4,0))
        tempOut[iTime,cellY,cellX] = temp[iTime,iHeight,iY,iX]



"""
This is an example program which writes some 4D pressure and
temperatures.
The companion program pres_temp_4D_rd.py shows how
to read the netCDF data file created by this program.

This example demonstrates the netCDF Python API.
It will work either with the Scientific Python NetCDF version 3 interface
(http://dirac.cnrs-orleans.fr/ScientificPython/)
of the 'classic' version of the netCDF4 interface.
(http://netcdf4-python.googlecode.com/svn/trunk/docs/netCDF4_classic-module.html)
To switch from one to another, just comment/uncomment the appropriate
import statements at the beginning of this file.

Jeff Whitaker <jeffrey.s.whitaker@xxxxxxxx> 20070202
"""

# the netCDF variable will be nrecs x nlevs x nlats x nlons.
nrecs = 2; nlevs = 2; nlats = 142; nlons = 311
#1*284*622

#lonOut = np.arange(-18.75, 59.00, 0.125)
#latOut = np.arange(24.75, 60.25, 0.125)

# open a new netCDF file for writing.
ncfile = Dataset('pres_temp_4D.nc','w')

# latitudes and longitudes of grid
#lats_out = 0.0 + 5.0*arange(nlats,dtype='float32')
#lons_out = 0.0 + 5.0*arange(nlons,dtype='float32')
lats_out = arange(start=24.75, stop=60.25, step=0.25,dtype='float32')
lons_out = arange(start=-18.75, stop=59.0, step=0.25,dtype='float32')

print len(lats_out)
print len(lons_out)

# output data.
press_out = 900. + arange(nlevs*nlats*nlons,dtype='float32') # 1d array
press_out.shape = (nlevs,nlats,nlons) # reshape to 2d array
#temp_out = 9. + arange(nlevs*nlats*nlons,dtype='float32') # 1d array
temp_out = tempOut
temp_out.shape = (nlevs,nlats,nlons) # reshape to 2d array

# create the lat and lon dimensions.
ncfile.createDimension('latitude',nlats)
ncfile.createDimension('longitude',nlons)

# create level dimension.
ncfile.createDimension('level',nlevs)

# create time dimension (record, or unlimited dimension)
ncfile.createDimension('time',None)

# Define the coordinate variables. They will hold the coordinate
# information, that is, the latitudes and longitudes.
# Coordinate variables only given for lat and lon.
lats = ncfile.createVariable('latitude',dtype('float32').char,('latitude',))

# Assign units attributes to coordinate var data. This attaches a
# text attribute to each of the coordinate variables, containing the
# units.
lats.units = 'degrees_north'
lons.units = 'degrees_east'

# write data to coordinate vars.
lats[:] = lats_out
lons[:] = lons_out

# create the pressure and temperature variables

# set the units attribute.
press.units =  'hPa'
temp.units = 'celsius'

# write data to variables along record (unlimited) dimension.
# same data is written for each record.
for nrec in range(nrecs):
   press[nrec,:,::] = press_out
   temp[nrec,:,::] = temp_out

# close the file.
ncfile.close()
print '*** SUCCESS writing example file pres_temp_4D.nc'

• Follow-Ups:
  • Re: [netcdfgroup] how to output a netcdf with an irregular grid
    • From: Signell, Richard