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

Fri, Jun 12, 2009 7:54 AM #

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

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]

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

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.

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
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.9 (MingW32)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iEYEARECAAYFAkoybCgACgkQzZ3g4YwFFgZ/+QCeOSkqOSQ7sqdCrLuGZtbkqZib
LAAAn3b5zqpAB1F/nGZDZw9HrhGdA95P
=HAAU
-----END PGP SIGNATURE-----

Edzer Pebesma

Fri, Jun 12, 2009 8:02 AM #

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

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.pebesma at wwu.de

Anne Ghisla

Fri, Jun 12, 2009 9:25 AM #

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Edzer Pebesma ha scritto:

hello all,

@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)[,
        2])
    convex_hull = c(convex_hull, convex_hull[1])
    d = Polygon(obj[convex_hull, ])
    new_data = spsample(d, 5000, type = "regular")
    gridded(new_data) = TRUE
    attr(new_data, "proj4string") <-obj at proj4string
    return(new_data)
}

Totally agree. The tests on projection are mandatory.

best regards,
Anne

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

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.9 (MingW32)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org

iEYEARECAAYFAkoygXYACgkQzZ3g4YwFFgZgDACeL7kQFjQGNqgz4mePX/nnCUkU
bzUAn37iRTmTsYAGiEWIOFWgM5P3UszT
=HqBt
-----END PGP SIGNATURE-----