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

[Anne Ghisla]

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

Paul Hiemstra ha scritto:
> 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.

Hi Paul, Edzer, list,

great news, thanks Paul for this improvement. autoKrige() use is now
more straightforward, I look forward to upgrade the package.
In the specific case of v.krige, anyway, I'll keep using a
self-generated grid from GRASS region and resolution.

> 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).

Useful information, I've missed it in previous threads.

best regards,
Anne

> 
> 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
>>>     
>>
>>   
> 
> 

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

iEYEARECAAYFAkpYiDYACgkQzZ3g4YwFFgY5rQCeKPqSz0T5VPCbeVEpPwbP0T48
TOYAn1v2cjFrnUHKOi6QlvhnXw3vOr93
=hVIx
-----END PGP SIGNATURE-----