[R-sig-Geo] [SoC v.autokrige2.py] autoKrige() and projected data

[Paul Hiemstra]

Hi Edzer, Anne and r-sig-geo,

I adapted the code of automap in reaction to this problem and included 
it in the new version of automap taht should be released somewhere next 
week.

from the documentation of autoKrige():

'autoKrige' performs some checks on the coordinate systems of
     'input_data' and 'new_data'. If one of both is 'NA', it is
     assigned the projection of the other. If they have different
     projections, an error is raised. If one of both has a
     non-projected system (i.e. latitude-longitude), an error is
     raised. This error is raised because 'gstat does use spherical
     distances when data are in geographical coordinates, however the
     usual variogram models are typically not non-negative definite on
     the sphere, and no appropriate models are available' (Edzer
     Pebesma on r-sig-geo).

This code illustrates the new functionality:

library(automap)
library(rgdal)
data(meuse)
coordinates(meuse) = ~x+y
data(meuse.grid)
gridded(meuse.grid) = ~x+y
proj4string(meuse) = CRS("+init=epsg:28992")

# new_data gets projection of input_data
autoKrige(log(zinc)~1, meuse, meuse.grid)

proj4string(meuse.grid) = CRS("+init=epsg:28992")
# Reproject to latlong
meuse_ll = spTransform(meuse, CRS("+init=epsg:4238"))
# Reproject to another projected system (LAEA)
meuse2 = spTransform(meuse, CRS("+init=epsg:3034"))
# Reproject to latlong
meuse.grid_ll = spTransform(meuse.grid, CRS("+init=epsg:4238"))

autoKrige(log(zinc)~1, meuse, meuse.grid)    # Should work
autoKrige(log(zinc)~1, meuse_ll, meuse.grid) # Should fail
autoKrige(log(zinc)~1, meuse, meuse.grid_ll) # Should fail
autoKrige(log(zinc)~1, meuse2, meuse.grid)   # Should fail

cheers,
Paul

Edzer Pebesma wrote:
> If newdata is created automatically (and this makes sense), it should
> recieve the projection information of data.
>
> On another issue: what happens if both are provided but have different
> projection attributes?
>
> And what happens if data has long/lat coordinates (is.projected(data) ==
> FALSE)? In that case I would prefer an error -- the current
> interpolation methods (covariances) used are not appropriate.
>
> Best regards,
> --
> Edzer
>
> Anne Ghisla wrote:
>   
>> Hello Paul, Edzer,
>>
>> the integration of automap package in v.autokrige2 is quite completed
>> [0], thanks Paul for this clean interface to gstat!
>> Sadly, I'm getting an error when autoKrige() receives no new_data
>> argument. Examples in documentation all involve meuse dataset, that has
>> no projection information. I browsed the Web for this situation with
>> little success. I don't know precisely where the error occurs, so I send
>> this report to Edzer too.
>>
>> Here is the code and a traceback of the error:
>>
>> ---code---
>> require(spgrass6)
>> require(automap)
>>
>> # rs in a vector map containing locations with relative elevation value,
>> in GRASS spearfish60 dataset. Obtained following this tutorial [1]
>>     
>>> rsmap <- readVECT6('rs', type='point')
>>>       
>> Exporting 300 points/lines...
>>  100%
>> 300 features written
>> OGR data source with driver: ESRI Shapefile
>> Source: "/home/anne/gis/spearfish60/PERMANENT/.tmp/galadriel", layer: "rs"
>> with  300  rows and  2  columns
>> Feature type: wkbPoint with 2 dimensions
>>
>>     
>>> str(rsmap)
>>>       
>> Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
>> [cut]
>>   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
>>   .. .. ..@ projargs: chr " +proj=utm +zone=13 +ellps=clrk66
>> +datum=NAD27 +units=m +no_defs
>> +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat"
>>
>> # not shown: meuse object has the same structure and contents but
>> projargs==NA.
>>
>>     
>>> autoKrige(elev~1, rsmap)
>>>       
>> Errore in predict.gstat(g, newdata = newdata, block = block, nsim =
>> nsim,  :
>>   var1 : data item in gstat object and newdata have different coordinate
>> reference systems
>>
>> Then I had a look into autoKrige code and found these lines:
>>
>>     if (missing(new_data))
>>         new_data = create_new_data(input_data)
>>     krige_result = krige(formula, input_data, new_data,
>> variogram_object$var_model,
>>         block = block, ...)
>>
>> My guess is that create_new_data() creates non-projected grid, and
>> krige() calls predict.gstat() with input_ and new_data, checks if they
>> have the same projection and in this case fails.
>>
>> In v.autokrige.py code, Mathieu Grelier workarounded this autogeneration
>> of new_data creating a correctly projected grid and passing it as
>> parameter.
>>
>> Shall I do the same or is it possible to get create_new_data() to keep
>> input_data's projection?
>>
>> many thanks in advance,
>> Anne
>>
>> [0] http://grass.osgeo.org/wiki/V.autokrige_GSoC_2009
>> [1] http://casoilresource.lawr.ucdavis.edu/drupal/node/438
>>     
>
>   


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul