[R-sig-Geo] Kriging and datum shifts

[Jeroen Steenbeek]

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
>
>