Dominik
Dominik

Reputation: 792

set Coordinate System when writing netcdf

I have a SpatialPointsDataFrame in R that I want to save as netCDF. I could use some help with formatting.

> str(swe)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 3487 obs. of  5 variables:
  .. ..$ site   : Factor w/ 6 levels "Dry Lake","Joe Wright",..: 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ depth  : num [1:3487] 151 157 138 155 145 ...
  .. ..$ density: num [1:3487] 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 ...
  .. ..$ swe.obs: num [1:3487] 0.56 0.582 0.512 0.572 0.535 ...
  .. ..$ date   : Date[1:3487], format: "2008-04-04" "2008-04-04" ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:3487, 1:2] -107 -107 -107 -107 -107 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "utm.e" "utm.n"
  ..@ bbox       : num [1:2, 1:2] -107.9 37.8 -104.1 43.8
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "utm.e" "utm.n"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr "+proj=longlat +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0"

I have a date field that is not regular so I think I can't used it as a dimension. These are field measurements of snow so the surveys are approximately once a month in two different years. I do however need to include the date and the site name for each measurement. I also need to preserve the projection information. Below is what I have so far:

require(ncdf4)
varlat=ncdim_def(name='latitude',units='deg',vals=coordinates(swe)[,2])
varlong=ncdim_def(name='longitude',units='deg',vals=coordinates(swe)[,1])
varswe=ncvar_def(name='swe.obs',units='meters',dim=list(varlat,varlong),missval=-9999)
varsite=ncvar_def(name='site',units='site',dim=list(varlat,varlong),missval=-9999)
vardate=ncvar_def(name='date',units='day',dim=list(varlat,varlong),missval=-9999)
new.nc=nc_create('snow.survey.nc',vars=list(varswe,varsite,vardate))

but when I try to populate values with:

 ncvar_put(new.nc,varid=varswe,vals=swe$swe.obs)
Error in ncvar_put(new.nc, varid = varswe, vals = swe$swe.obs) : 
  ncvar_put: error: you asked to write 12159169 values, but the passed data array only has 3487 entries!

here is new.nc:

> new.nc
[1] "File snow.survey.nc (NC_FORMAT_CLASSIC):"
[1] ""
[1] "     3 variables:"
[1] "        float swe.obs[latitude,longitude]   "
[1] "            units: meters"
[1] "            _FillValue: -9999"
[1] "        float site[latitude,longitude]   "
[1] "            units: site"
[1] "            _FillValue: -9999"
[1] "        float date[latitude,longitude]   "
[1] "            units: day"
[1] "            _FillValue: -9999"
[1] ""
[1] "     2 dimensions:"
[1] "        latitude  Size:3487"
[1] "            units: deg"
[1] "            long_name: latitude"
[1] "        longitude  Size:3487"
[1] "            units: deg"
[1] "            long_name: longitude"

Also, How do I define the projection string in the netcdf file? in this case '+proj=longlat +datum=NAD83' so another user knows what they are dealing with? Thanks

UPDATE1: I tried this based on @Spacedman but I still get an error. Do I need to make a regular grid of the extents of my coordinates and then try to use that as the extents of my dimension? I'd have ot fill in NA everywhere I don't have a data point. I also don't have evenly spaced measurements.

coords=cbind(coordinates(swe)[,1],coordinates(swe)[,2])
nccoords=ncdim_def(name='coords',units='site',vals=coords)
varswe=ncvar_def(name='swe.obs',units='meters',dim=nccoords,missval=-9999,longname='survey points')

create new 'new.nc' simplified to one variable at the moment.

new.nc=nc_create('snow.survey.nc',vars=varswe)
ncvar_put(new.nc,varid=varswe,vals=swe$swe.obs)
Error in ncvar_put(new.nc, varid = varswe, vals = swe$swe.obs) : 
  ncvar_put: error: you asked to write 6974 values, but the passed data array only has 3487 entries!

Upvotes: 2

Views: 1165

Answers (1)

Spacedman
Spacedman

Reputation: 94202

You've got your dimensions in a twist. Your lat-long data should be saved in a 3487x2 dimension variable, and your attributes in 3487x1 variables. You've created things (varswe, varsite...) that are 3487x3487 in dimension. This is not a grid.

If your points are on a grid (it would have to be a 317x11 grid if complete) then you need to get the unique lat and long values and create dimensions of length 317 and 11 for those.

Upvotes: 2

Related Questions