Kriging and datum shifts
Thank you Mike,
I have indeed all the info that I need in the SPDF: lat, lon, CRS. The
irregular distribution of data is not projection-related, but is
prohibitive to using simple rasterization operations because source data
and target cells do not always overlap. Hence the need for kriging (or
some other form or interpolation). Rasterization would have made life so
much easier... I will find a way, your suggestions are helpful.
As a though: I am just surprised to find that the sp package does not
seem to notice that two global datasets with different longitude offsets
entirely overlap. In my case, kriging the point dataset from -280 to +80
lon to a raster from -180 to 180 lon leaves the raster data +80 and
onward blank, while in fact the point dataset has provided data for the
entire raster range - just at a different lon offset. It is an
understandable but unfortunate shortcoming shared by many (GIS) toolkits.
I can work around this by modifying the longitude of the SPDF points:
lon <- (lon + 180) %% 360 - 180
That is valid but rather hack.... Helmet time.
This leads me to a second issue: I suspect that when operating idw on a
global dataset, interpolation is also not properly executed around the
wrapping longitude boundary (e.g., I suspect that kriging is not applied
to the boundary line at +180 for a raster that spans -180 to 180
longitude, creating noticeable artifacts)
Best,
Jeroen
On 2013-05-27 6:03 PM, Michael Sumner wrote:
For a simple approach I would just do on the raw coordinates in the
SPDF. You could operate on the NetCDF file itself by reading from it
directly as a matrix with vectors, but since you mention the SPDF why
not just do this:
x <- as.data.frame(spdf)
## figure out what "longitude" and "latitude" are called here
lon <- x$longitude
x$longitude <- ifelse(lon < -180, lon + (-180 - -280), lon)
## reconstruct, also with what "longitude" and "latitude" are named
xr <- SpatialPointsDataFrame(x[, c("longitude", "latitude")], x,
proj4string = CRS(proj4string(spdf)))
(Though I doubt you have a valid CRS string in your original object
already, the above will work. To give it one replace the proj4string
call with "+proj=longlat +datum=WGS84").
Now proceed with xr instead.
There is rotate() in the raster package, though it's only for one
special case. You can easily crop and re-merge the two parts as
rasters with that package too, and you can do the same with and
subsetting and recombining in sp, but in this case I would go as far
back to the source first. R has excellent tools for dealing with data
that need some kind of pre-preparation before they are ready for these
spatial tools. As an example, sp will baulk at this
SpatialPointsDataFrame(cbind(c(-280, -100), c(-40, -60)),
data.frame(x = 1:2), proj4string = CRS("+proj=longlat"))
Since you are reading data from NetCDF and they are on an irregular
grid I would also spend time becoming familiar with that by reading
from it with generic tools like ncdf or ncdf4 or RNetCDF, R's image()
list structure will handle irregular coordinates in one or two
dimensions and sometimes that can save you from confusion or pain.
Also, sometimes NetCDF data are using some underlying map projection
and if you can discover that rather than work in the irregular lon/lat
you'll be able to cast more directly to these rigid grid/raster
objects. Mercator is commonly used for global grids and provides
irregular latitudes for example. Can you point to one of these files?
Cheers, Mike.
On Tue, May 28, 2013 at 12:50 AM, Jeroen Steenbeek <drmbongo at gmail.com> wrote:
Dear list, I am converting NetCDF data, placed on a global grid with irregular latitude values, to a regular WGS84 grid using IDW. It works, but I have one remaining issue for which I can use advice. The original data is oriented -280, +80 longitude, and it needs to be interpolated to -180, 180 longitude. How do I do this with the Spatial* classes? I can do this either before interpolation on the original SpatialPointsDataFrame that is read from the NetCDF files, or I can perform this shift on the SpatialGridDataFrame resulting from the IDW operation. But how to do this in R? spCBind sounds like a good candidate, but will this correctly handle the georeferencing? Thanks in advance, Jeroen
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo