[R-sig-Geo] how to write a netCDF with irregular grid

[Michael Nolde]

Hello there,

I've got a netCDF file which I want to read, modify, and save, which works
well with 'ncdf4'-package.
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"', which,
as far as I understood, make my viewer application (Panoply) make the data
appear in the correct locations.
However, I am not able to hand this argument to the 'ncvar_def' method (the
error is: 'unused argument'), so at the moment I use only one dimension of
the 'lat' and 'lon' variables and hand them to the 'ncvar_def'-method, which
produces an distorted (regular) output.

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 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 source code so far:
---------------------------------

library(ncdf4)
library(maptools)
library(sp)

ncdfname <- "TMAX_2M.nc"
ncdffile <- nc_open(ncdfname)
latorig <- ncvar_get(ncdffile,"lat",verbose=F)
lonorig <- ncvar_get(ncdffile,"lon")
timeorig <- ncvar_get(ncdffile,"time",verbose=F)
tunits <- ncatt_get(ncdffile,"time","units")

# get global attributes
titleorig <- ncatt_get(ncdffile,0,"title")
institutionorig <- ncatt_get(ncdffile,0,"institution")
datasourceorig <- ncatt_get(ncdffile,0,"source")
referencesorig <- ncatt_get(ncdffile,0,"references")
historyorig <- ncatt_get(ncdffile,0,"history")
Conventionsorig <- ncatt_get(ncdffile,0,"Conventions")

temp.array <- ncvar_get(ncdffile,"TMAX_2M")
dlname <- ncatt_get(ncdffile,"TMAX_2M","long_name")
dunits <- ncatt_get(ncdffile,"TMAX_2M","units")
fillvalue <- ncatt_get(ncdffile,"TMAX_2M","_FillValue")
temp.array[temp.array==fillvalue$value] <- NA
temp.array[temp.array==NULL] <- NA
nc_close(ncdffile)
temp.slice <- temp.array[,,1]


lonorig1d = lonorig[,1]
latorig1d = latorig[1,]

lonorig1d_sort = sort(lonorig1d, decreasing=FALSE)
latorig1d_sort = sort(latorig1d, decreasing=FALSE)

dim_lon  <- ncdim_def('x', 'degrees_east', as.double(lonorig1d_sort),
create_dimvar=TRUE)
dim_lat  <- ncdim_def('y', 'degrees_north', as.double(latorig1d_sort),
create_dimvar=TRUE)

dim_height  <- ncdim_def('height_2m', 'height', c(0:0))
dim_time <- ncdim_def('time', "seconds since 1979-01-01 00:00:00",
calendar="proleptic_gregorian", unlim=T, c(364:0))

var_out_lon <- ncvar_def('lon', 'degrees_east', list(dim_lon),
longname="longitude")
var_out_lat <- ncvar_def('lat', 'degrees_north', list(dim_lat),
longname="latitude")
var_out_timebnds <- ncvar_def('time_bnds', "seconds since 1979-01-01
00:00:00", list(dim_time))

var_out <- ncvar_def('TMAX_2M', 'K', list(dim_lon,dim_lat,dim_time), 9.e20,
longname="2m maximum temperature")

ncid_out <- nc_create('updated_netcdf.nc', list(var_out, var_out_lon,
var_out_lat))

ncvar_put(ncid_out, var_out, temp.array[,,])
nc_close(ncid_out)