[R-sig-Geo] [SoC v.autokrige2.py] autoKrige() and projected data
edzer.pebesma at uni-muenster.de
Fri Jun 12 17:02:55 CEST 2009
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
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.
Anne Ghisla wrote:
> Hello Paul, Edzer,
> the integration of automap package in v.autokrige2 is quite completed
> , 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:
> # rs in a vector map containing locations with relative elevation value,
> in GRASS spearfish60 dataset. Obtained following this tutorial 
> > rsmap <- readVECT6('rs', type='point')
> Exporting 300 points/lines...
> 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
> ..@ 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
> > 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,
> 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
> Shall I do the same or is it possible to get create_new_data() to keep
> input_data's projection?
> many thanks in advance,
>  http://grass.osgeo.org/wiki/V.autokrige_GSoC_2009
>  http://casoilresource.lawr.ucdavis.edu/drupal/node/438
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de
More information about the R-sig-Geo