Interpolation

Interpolation

As an example we will interpolate the altitude of the peak shapefile.

library(rgdal)
srtm <- readGDAL("./results/srtm.tif")

## ./results/srtm.tif has GDAL driver GTiff 
## and has 481 rows and 289 columns

str(srtm@proj4string)

## Formal class 'CRS' [package "sp"] with 1 slot
##   ..@ projargs: chr "+proj=utm +zone=44 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

peaks <- readOGR(dsn = "./data/DSL250-Shp/", layer = "Peaks")

## OGR data source with driver: ESRI Shapefile 
## Source: "./data/DSL250-Shp/", layer: "Peaks"
## with 21 features
## It has 6 fields

str(peaks)

## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  21 obs. of  6 variables:
##   .. ..$ AREA     : num [1:21] 0 0 0 0 0 0 0 0 0 0 ...
##   .. ..$ PERIMETER: num [1:21] 0 0 0 0 0 0 0 0 0 0 ...
##   .. ..$ PEAKS_   : num [1:21] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..$ PEAKS_ID : num [1:21] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..$ NAME     : Factor w/ 20 levels "ALAGALLA","BIBLE ROCK",..: 15 19 8 6 NA 1 16 7 2 5 ...
##   .. ..$ ALTITUDE : int [1:21] 766 254 293 534 525 1033 2243 661 798 1358 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:21, 1:2] 186782 182254 213122 240116 149359 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
##   ..@ bbox       : num [1:2, 1:2] 148456 126653 285930 373827
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=tmerc +lat_0=7.000480277777778 +lon_0=80.77171111111112 +k=0.9999238418 +x_0=200000 +y_0=200000 +a=6377276.345 +b=6356075"| __truncated__

peaks <- spTransform(x = peaks, CRSobj = srtm@proj4string)

We will start with a determinstic interpolation approach that you should be familiar with: Inverse Distance Weighting (IDW).

library(gstat)
peak_idw <- idw(peaks@data$ALTITUDE ~ 1, peaks, srtm, maxdist = 200000, idp = 2)

## [inverse distance weighted interpolation]

#spplot(peak_idw)
plot(raster(peak_idw))
#contour(peak_idw, add=TRUE)
points(peaks,pch=20,cex=.4)

Now, let us try a statistical interpolation approach, i.e. simple Kriging.

library(gstat)

plot(variogram(peaks@data$ALTITUDE ~ 1, peaks, cloud = TRUE))

plot(variogram(peaks@data$ALTITUDE ~ 1, peaks))

plot(variogram(peaks@data$ALTITUDE ~ 1, peaks, alpha = c(0,45,90,135)))

vt <- variogram(peaks@data$ALTITUDE ~ 1, peaks)

show.vgms()

v.fit <- fit.variogram(vt,
    vgm(nugget = 0, model = "Gau", psill = 800000, range = 40000))
v.fit2 <- fit.variogram(vt,
    vgm(nugget = 0, model = "Sph", psill = 800000, range = 40000))
plot(vt,v.fit)

plot(vt,v.fit2)

peak_k <- krige(peaks@data$ALTITUDE ~ 1, peaks, srtm, v.fit,nmin = 3, maxdist = 300000, nmax = 7)

## [using ordinary kriging]

peak_k2 <- krige(peaks@data$ALTITUDE ~ 1, peaks, srtm, v.fit2,nmin = 3, maxdist = 300000, nmax = 7)

## [using ordinary kriging]

par(mfrow = c(1,2))
plot(raster(peak_k),main="Kriging using Gaussian fit")
points(peaks,pch=20,cex=.4)
plot(raster(peak_k2),main="Kriging using Spherical fit")
points(peaks,pch=20,cex=.4)