[SoC v.autokrige2.py] autoKrige() and projected data

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