Since your question (and possible answers) involve code, I suggest
submitting this question to stackoverflow.com or gis.stackexchange.com
and tagging with "netcdf" and "python". It typically works a lot
better than e-mail for these type of questions.
On Thu, Feb 27, 2014 at 1:18 PM, Michael Nolde
<nolde@xxxxxxxxxxxxxxxxxxxxxx> wrote:
> Dear list members,
>
> I have been trying to write a netcdf with an irregular grid for weeks now,
> but failed. I used 'Python' with 'python-netcdf4' package, as well as 'R'
> with
> 'ncdf4'-package. I also installed 'cdo' with netcdf-support and tried to
> reproject
> 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
> values from the irregular grid to it, which costs resolution and takes
> considerably
> more disk space (source code attached below). Of course, I searched the web
> up
> and down, but could not find anything related to that problem. Could anyone
> please help me?
> Any way to create a netcdf file with an irregular grid would be highly
> appreciated.
>
>
>
> 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.
> Below I'm posting the structure of my input file, which I hope to reproduce,
> 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',))
> lons =
> ncfile.createVariable('longitude',dtype('float32').char,('longitude',))
>
> # 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
> press =
> ncfile.createVariable('pressure',dtype('float32').char,('time','level','latitude','longitude'))
> temp =
> ncfile.createVariable('temperature',dtype('float32').char,('time','level','latitude','longitude'))
>
> # 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'
>
> _______________________________________________
> netcdfgroup mailing list
> netcdfgroup@xxxxxxxxxxxxxxxx
> For list information or to unsubscribe, visit:
> http://www.unidata.ucar.edu/mailing_lists/
--
Dr. Richard P. Signell (508) 457-2229
USGS, 384 Woods Hole Rd.
Woods Hole, MA 02543-1598
