Estimated soil moisture

In a previous document, we estimated soil moisture from Landsat-based multispectral indices and SRTM-based topographic variables. In order to conduct accuracy assessment of such an estimation, we will use here SMAP data as reference. First of all, let’s read the estimated soil moisture. Note that this dataset is 30 m spatial resolution.

library(raster)
## Loading required package: sp
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
sm <- raster("C:\\Users\\Iván Lizarazo\\Documents\\UN\\datos\\soilmoisture\\piedemonte\\rf_sm2015_06_10.tif")

sm
## class       : RasterLayer 
## dimensions  : 2545, 3055, 7774975  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1028016, 1119666, 934053, 1010403  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : C:\Users\Iván Lizarazo\Documents\UN\datos\soilmoisture\piedemonte\rf_sm2015_06_10.tif 
## names       : rf_sm2015_06_10 
## values      : 27.54736, 43.82128  (min, max)
plot(sm)

SMAP data for 19 June 2015

SMAP data can be read as a raster object using the raster package. Note that this is a Level 3 dataset with 3 km spatial resolution.

smap <- raster("C:\\Users\\Iván Lizarazo\\Documents\\UN\\datos\\soilmoisture\\piedemonte\\test\\smap_19Jun2015.tif")

smap <- smap * 100

plot(smap)

W <- shapefile("C:\\Users\\Iván Lizarazo\\Documents\\UN\\datos\\soilmoisture\\piedemonte\\WGuatiquia_Negro.shp")

plot(W, add=TRUE)

It’s obvious that a comparison cannot be made without SMAP interpolation. We will use here a IDW technique for conducting such a task.

Converting SMAP to spatial points dataframe

psmap <- as(smap, 'SpatialPointsDataFrame')

psmap
## class       : SpatialPointsDataFrame 
## features    : 154 
## extent      : -73.8509, -72.9169, 3.819202, 4.737002  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       :   smap_19Jun2015 
## min values  : 11.1691551917669 
## max values  :  52.808390835191

Visualizing data

library(tmap)

library(leaflet)

tm_shape(W) + tm_borders(col="green") +
  tm_shape(psmap) + 
  tm_dots(col="smap_19Jun2015", palette="RdBu",
        auto.palette.mapping =FALSE,
             title="SMAP soil moisture (in %)",
  size=0.1) +

  tm_legend(legend.outside=TRUE)

IDW interpolation of SMAP dataset

The IDW output is a raster. This requires that we first create an empty raster grid, then interpolate the precipitation values to each unsampled grid cell. An IDW power value of 2 (idp=2.0) will be used.

library(gstat) # Use gstat's idw routine
library(sp)    # Used for the spsample function

# Create an empty grid where n is the total number of cells
grd              <- as.data.frame(spsample(psmap, "regular", n=5000000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Create SpatialPixel object
fullgrid(grd)    <- TRUE  # Create SpatialGrid object

# Add P's projection information to the empty grid
proj4string(grd) <- proj4string(psmap)

# Interpolate the grid cells using a power value of 2 (idp=2.0)
P.idw <- gstat::idw(smap_19Jun2015 ~ 1, psmap, newdata=grd, idp=2.0)
## [inverse distance weighted interpolation]
# Convert to raster object then clip to ROI
r       <- raster(P.idw)
#r.m     <- mask(r, W)

# Plot
tm_shape(W) + tm_borders(col="green") +
tm_shape(r) + 
  tm_raster(n=10,palette = "RdBu", auto.palette.mapping = FALSE,
            title="Interpolated soil moisture \n(in inches)") + 
  tm_shape(psmap) + tm_dots(size=0.2) +
  tm_legend(legend.outside=TRUE)

Now, we need to put the r raster into the sm raster’s coordinate reference system.

rproj <- projectRaster(r,sm)

rf <- mask(rproj, sm)
### check both raster are similar
rf
## class       : RasterLayer 
## dimensions  : 2545, 3055, 7774975  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1028016, 1119666, 934053, 1010403  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : var1.pred 
## values      : 13.67086, 49.35837  (min, max)
sm
## class       : RasterLayer 
## dimensions  : 2545, 3055, 7774975  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1028016, 1119666, 934053, 1010403  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : C:\Users\Iván Lizarazo\Documents\UN\datos\soilmoisture\piedemonte\rf_sm2015_06_10.tif 
## names       : rf_sm2015_06_10 
## values      : 27.54736, 43.82128  (min, max)

Calculating RMSD and MAD between estimated soil moisture and SMAP data

### write RMSE and R2 functions
calculate.RMSE <- function (x,y){
  # x is SMAP data, y is the predicted soil moisture
  minus<-x-y
  square<-minus^2
  mean<-cellStats(square, 'mean')
  RMSE <- sqrt(mean)
  return(RMSE)
}
#
# Function that returns Mean Absolute Error
calculate.MAE <- function(x,y)
{
  minus<-x-y
  absolute <- abs(minus)
  media <- cellStats(absolute, 'mean')
  MAE <- media
  return(MAE)
}

### use the functions
RMSEpredict <- calculate.RMSE(rf, sm)
paste("The RMSD is", format(round(RMSEpredict,2), nsmall = 2))
## [1] "The RMSD is 6.15"
MAEpredict <- calculate.MAE(rf, sm)
paste("The MAD is", format(round(MAEpredict,2), nsmall = 2))
## [1] "The MAD is 4.78"

Soil moisture visual comparison

library(leaflet)

           
pal <- colorNumeric(palette="RdYlBu", domain = c(0, 100),
 na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(rf, colors = pal, opacity = 1.0) %>%
  addLegend(pal = pal, values = values(rf),
    title = "SMAP interpolated soil moisture - June 2015")
library(leaflet)

#pal <- colorNumeric(palette="RdYlBu", domain = c(0, 100), n=10,
 # na.color = "transparent")

leaflet() %>% addTiles() %>%
  addRasterImage(sm, colors = pal, opacity = 1.0) %>%
  addLegend(pal = pal, values = values(sm),
    title = "Estimated soil moisture - June 2015")