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 > +nadgri...@conus,@alaska,@ntv2_0.gsb,@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 -- Edzer Pebesma 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.pebe...@wwu.de _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo