Importation

Read dataset from gstat

data(sic97)

mydata<-st_as_sf(sic_full)
mydataobs<-st_as_sf(sic_obs)
mydata$obs<- mydata$ID %in% mydataobs$ID

# From internet (data projection)
crs.suisse<-"+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1  +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +units=m +no_defs +x_0=-55000 +y_0=14000"
st_crs(mydata)<-crs.suisse
myobs<-mydata[mydata$obs,]

Download elevation

if (!("SRTMsuisse.tif" %in% list.files(file.path(getwd(),"datasets"))))
{
  elev<-NULL
  for (LONG in c("006","007","008","009","010"))
    for (LAT in c("45","46","47"))
    {
      TILE<-paste0("N",LAT,"E",LONG)
      urlzip<-paste0("http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/",TILE,".hgt.zip")
      download.file(url=urlzip,destfile=paste0(TILE,".hgt.zip"),mode="wb")
      unzip(zipfile=paste0(TILE,".hgt.zip"))
      srtm<-raster(paste0(TILE,".hgt"))
      if (is.null(elev)) elev<-srtm
      if (!is.null(elev)) elev<-merge(elev,srtm)
    }
  
  # reproject elev to swiss CRS
  elev<-projectRaster(elev, crs=crs.suisse)
}
if ("SRTMsuisse.tif" %in% list.files(file.path(getwd(),"datasets"))) 
  elev<-raster(file.path(getwd(),"datasets","SRTMsuisse.tif")) # 

# Check that elev is a raster layer
elev
## class      : RasterLayer 
## dimensions : 3678, 6177, 22719006  (nrow, ncol, ncell)
## resolution : 64, 92.6  (x, y)
## extent     : -168975.9, 226352.1, -203412.9, 137169.9  (xmin, xmax, ymin, ymax)
## crs        : +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=-55000 +y_0=14000 +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +units=m +no_defs 
## source     : C:/Users/mistea/Documents/OpenSpat/Gembloux juin 2018/Practical session co-kriging/datasets/SRTMsuisse.tif 
## names      : SRTMsuisse 
## values     : -7.423397, 4672  (min, max)
# confirm that the DEM data covers all our precipitation dataset:
par(mfrow=c(1,1),mar=rep(0,4))
plot(elev)
mydata %>% st_coordinates() %>% points()

raster::extract elevation values from DEM at our locations

e<-extract(elev,st_coordinates(mydata))
mydata$elev<-e

Check elevation and compare with Topography

# not necessary: just to confirm that the elevation values make sense:
# plot elevation values over Topography
mapviewOptions(basemaps="OpenTopoMap")
mapview(mydata,zcol="elev",lwd=0)

NOTE: there are NA values in elevation: possibly because the SRTM has missing values at some locations

sum(is.na(mydata$elev)) # 1 NA values
## [1] 1
# to remove those locations from the data set:
mydata<-mydata[!is.na(mydata$elev), ]
myobs<-mydata[mydata$obs,]

Building a dataset for experimental Variograms - Cross Variograms

suisse.gs <- gstat(NULL,id="rainfall", formula = rainfall~1,  data=mydata) # set = list(nocheck = 1),
suisse.gs <- gstat(suisse.gs,id="elev", formula = elev~1, data=mydata) # set = list(nocheck = 1),

Model for Rainfall Variograms - Cross Variograms.

#Variograms
suisse.vg <- gstat::variogram(suisse.gs)
plot(suisse.vg, main='Rainfall - Elevation Variograms')

# Guess at a variogram model for each
# add the model to the gstat object
suisse.gs <- gstat(suisse.gs, model = vgm("Sph", nugget = 200, range = 90000, psill=15000), fill.all=TRUE) # set = list(nocheck = 1)

Automatic fit

# Cross-Variograms
suisse.fit <- gstat::fit.lmc(suisse.vg, suisse.gs, fit.lmc=TRUE) 
plot(suisse.vg, model=suisse.fit, 
     main="Fitted Variogram Models - Raw Data")

print(suisse.fit)
## data:
## rainfall : formula = rainfall`~`1 ; data dim = 466 x 4
## elev : formula = elev`~`1 ; data dim = 466 x 4
## variograms:
##                  model      psill range
## rainfall[1]        Nug    296.888     0
## rainfall[2]        Sph  14470.072 90000
## elev[1]            Nug  20831.476     0
## elev[2]            Sph 233133.770 90000
## rainfall.elev[1]   Nug   2253.722     0
## rainfall.elev[2]   Sph  -8021.999 90000
# to modify
# suisse.fit$model$elev$psill[2] <- 600000      # nugget for elevation
# suisse.fit$model$elev$range[2] <- 100000       # psill for elevation
#etc.

Prediction

suisse_nobs <- mydata[!mydata$obs,] # sf

# predict to the non observed points
cok <- predict(suisse.fit, newdata=suisse_nobs) 
## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]
# summarize predictions and their errors
summary(cok$rainfall.pred); summary(cok$rainfall.var)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   100.0   162.0   185.4   263.5   517.0
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -2.001e-11 -5.457e-12  0.000e+00 -2.528e-13  3.638e-12  2.547e-11

Interpolation Map

First, we need a grid.

chull<-suisse_nobs %>% st_union() %>% st_convex_hull() 
suisse.points<-suisse_nobs %>% st_make_grid(n=100) %>% st_intersection(chull) %>% st_centroid() %>% st_sf() #sf
suisse.grid<-suisse_nobs %>% st_make_grid(n=100)  %>% st_intersection(chull)

Estimation at the grid nodes

cok.grid <- predict(suisse.fit, newdata=suisse.points)
## Linear Model of Coregionalization found. Good.
## [using ordinary cokriging]
st_geometry(cok.grid)<-st_geometry(suisse.grid) # to assign polygons

Interpolation Map

mapviewOptions(basemaps="CartoDB.Positron")
mapview(cok.grid, zcol="rainfall.pred",lwd=0)+mapview(suisse_nobs,zcol="rainfall")

Kriging Variance

mapviewOptions(basemaps="CartoDB.Positron")
mapview(cok.grid, zcol="rainfall.var",lwd=0)+mapview(suisse_nobs,col.regions="red",cex=3)