[R-sig-Geo] [SoC v.autokrige2.py] autoKrige() and projected data
a.ghisla at gmail.com
Fri Jun 12 18:25:26 CEST 2009
-----BEGIN PGP SIGNED MESSAGE-----
Edzer Pebesma ha scritto:
> If newdata is created automatically (and this makes sense), it should
> recieve the projection information of data.
@Paul: I added the projection of obj to new_data before return. Please
test if it works as part of the package.
create_new_data <- function (obj)
convex_hull = chull(coordinates(obj)[, 1], coordinates(obj)[,
convex_hull = c(convex_hull, convex_hull)
d = Polygon(obj[convex_hull, ])
new_data = spsample(d, 5000, type = "regular")
gridded(new_data) = TRUE
attr(new_data, "proj4string") <-obj at proj4string
> 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.
Totally agree. The tests on projection are mandatory.
> Best regards,
> 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
>> 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
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.9 (MingW32)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org
-----END PGP SIGNATURE-----
More information about the R-sig-Geo