[R-sig-Geo] Kriging and datum shifts

[Edzer Pebesma]

Jeroen, sp does not have anything to do with kriging, so your comment is
confusing.

sp does create spatial objects from data.frame's, and may do some sanity
checks on coordinates in case they are longlat -- you have to specify
they are longlat, though.

If they are specified as longlat, and e.g. gstat::idw is used for
interpolation, it should use great circle distances as opposed to
Euclidian distances to figure out interpolation weights. If this is not
the case, I'd like to look at a simple, reproducible example. Kriging on
spherical data is done at your own risk, there are no checks that the
variogram model can be applied sanely on a sphere, as reported before.

On 05/28/2013 03:47 PM, Jeroen Steenbeek wrote:
> 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
>>
>>
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
83 33081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de