[R-sig-Geo] correct use of autoKrige and gstat with Digital Elevation Models

Stefano Sofia @tef@no@@ofi@ @ending from regione@m@rche@it
Thu Aug 16 09:21:15 CEST 2018


Dear R-sig-geo list users,
I am dealing with rainfall isohyet maps, I am using a DEM in order to get more accurate values.
In particular I am trying to compare Kriging method (through autoKrige) with the Inverse Distance method (through gstat).
Kriging code works fine, I think that I am not able to make a correct use of gstat.

In both cases,
- I load my rainfall cumulates in a data frame called df_prec:

Codice_sensore Cumulata Long_Cent Lat_Cent Quota
3134 177 13.22888 42.99361  1352
3135 98 13.19046 42.85843  1496
3136 40 12.81444 43.23840  1005
3137 201 13.30681 42.84369  1036
3138 145 13.32638 42.90215   980
3139 45 13.11180 42.79578  1575
3140 69 13.18538 42.91439  1800
3234 129 13.44506 42.74417   960

- I load the shape file and the DEM:

## LOAD SHAPEFILE
Shape_area5 <- readOGR(dsn="area5.shp", layer="area5")
Shape_area5_projection <- "+init=epsg:32633 +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
Shape_area5 <- spTransform(Shape_area5, CRS(Shape_area5_projection))
## LOADING DEM
ita_DEM <- getData('alt', country='ITA', mask=FALSE)
ita_DEM <- projectRaster(ita_DEM, crs = CRS("+init=epsg:32633"))

- I put the coordinates in the correct way, I evaluate the extracted elevation values and I create a new grid with these new values:

##
coordinates(df_prec) <- ~Long_Cent+Lat_Cent
proj4string(df_prec) <- CRS("+init=epsg:4326")
df_prec <- spTransform(df_prec, CRS("+init=epsg:32633"))
## EXTRACT THE ELEVATION VALUES TO MY POINTS; extract EXPECTS A RASTER INPUT
df_prec$ExtractedElevationValues <- extract(x=ita_DEM, y=df_prec)
## CREATE A NEW GRID
newgrid <- as(ita_DEM, "SpatialGridDataFrame")
names(newgrid) <- "ExtractedElevationValues"

- when I apply Kriging, everything works fine
## KRIGING
RegressionKriging_rainfall <- autoKrige(Cumulata~ExtractedElevationValues, df_prec, newgrid)
prec_rast <- raster(RegressionKriging_rainfall$krige_output)

- if I use gstat in the standard way

Inverse_Distance_rainfall <- gstat(formula=Cumulata~1, data=df_prec, locations=newgrid)

everything goes fine, but this code obviously does not take into consideration the elevation points. When I try to use the inverse distance method like Kriging:
Inverse_Distance_rainfall <- gstat(formula=Cumulata~ExtractedElevationValues, df_prec, newgrid)

I get the following error:

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function �proj4string� for signature �"NULL"�
  Calls: my_inversedistanceweighted_DEM -> gstat -> proj4string -> <Anonymous>


Where is my mistake?
Could please somebody help me and give me a hint for a correst use of gstat?
Thank you for your attention and your help
Stefano


         (oo)
--oOO--( )--OOo----------------
Stefano Sofia PhD
Area Meteorologica e  Area nivologica - Centro Funzionale
Servizio Protezione Civile - Regione Marche
Via del Colle Ameno 5
60126 Torrette di Ancona, Ancona
Uff: 071 806 7743
E-mail: stefano.sofia using regione.marche.it
---Oo---------oO----------------

________________________________

AVVISO IMPORTANTE: Questo messaggio di posta elettronica pu� contenere informazioni confidenziali, pertanto � destinato solo a persone autorizzate alla ricezione. I messaggi di posta elettronica per i client di Regione Marche possono contenere informazioni confidenziali e con privilegi legali. Se non si � il destinatario specificato, non leggere, copiare, inoltrare o archiviare questo messaggio. Se si � ricevuto questo messaggio per errore, inoltrarlo al mittente ed eliminarlo completamente dal sistema del proprio computer. Ai sensi dell�art. 6 della DGR n. 1394/2008 si segnala che, in caso di necessit� ed urgenza, la risposta al presente messaggio di posta elettronica pu� essere visionata da persone estranee al destinatario.
IMPORTANT NOTICE: This e-mail message is intended to be received only by persons entitled to receive the confidential information it may contain. E-mail messages to clients of Regione Marche may contain information that is confidential and legally privileged. Please do not read, copy, forward, or store this message unless you are an intended recipient of it. If you have received this message in error, please forward it to the sender and delete it completely from your computer system.

--
Questo messaggio  stato analizzato da Libra ESVA ed  risultato non infetto.
This message was scanned by Libra ESVA and is believed to be clean.


	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list