In villandre/ISMRAanalyses: Analyses for IS-MRA study
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
We start by loading the required libraries.
if ("package:ISMRAanalyses" %in% search()) {
detach("package:ISMRAanalyses", unload = TRUE)
}
library(ISMRAanalyses)
library(sp)
library(raster)
library(MODIS)
library(rgdal)
library(rgeos)
library(spacetime)
library(MRAinla)
library(RhpcBLASctl)
library(geoR)
library(maptools)
library(doParallel)
library(mapmisc)
library(xtable)
# Input the correct working directory for your system
# **The working directory should have two subfolders: "data" for datasets and "outputFiles" for model fitting outputs and graphs.**
setwd("/home/luc/INLAMRAfiles/INLAMRApaper1/realData")
RandomFields::RFoptions(cores = 1)
blas_set_num_threads(1)
omp_set_num_threads(1)
Now, we load all datasets required for producing the results presented in the article.
# Input in variable rawDataFilesLocation the folder name where MODIS data files are found
rawDataFilesLocation <- "/store/luc/rawDataFiles"
trainingDataValidationName <- load("data/mainDataCompleteMap_May18_24.Rdata")
trainingDataValidation <- get(trainingDataValidationName)
rm(list = trainingDataValidationName)
testDataValidationName <- load("data/testDataMay21_May18_24.Rdata")
testDataValidation <- get(testDataValidationName)
rm(list = testDataValidationName)
indiaAnalysisValidationName <- load("outputFiles/INLAMRA_validationAnalysis.Rdata")
indiaAnalysisValidation <- get(indiaAnalysisValidationName)
if (indiaAnalysisValidationName != "indiaAnalysisValidation") {
rm(list = indiaAnalysisValidationName)
}
SPDEresultName <- load("outputFiles/inlaFitForValidationAnalysis.Rdata")
SPDEresult <- get(SPDEresultName)
rm(list = SPDEresultName)
predictionDataMainName <- load("data/testDataMay28_May25_31_Larger.Rdata")
predictionDataMain <- get(predictionDataMainName)
rm(list = predictionDataMainName)
trainingDataMainName <- load("data/mainDataCompleteMap_May25_31_Larger.Rdata")
trainingDataMain <- get(trainingDataMainName)
rm(list = trainingDataMainName)
indiaAnalysisMainName <- load("outputFiles/indiaAnalysisMay28evenLarger_May25_31_VersionApril13.Rdata")
indiaAnalysisMain <- get(indiaAnalysisMainName)
rm(list = indiaAnalysisMainName)
We now produce the basic raster plots for the observed land surface temperatures, with cities and geopolitical boundaries overlaid.
# Naming convention: nnnnnnn.Ayyyyddd.h00v00.vvv.yyyydddhhmmss.
# nnnnnnn: Product name
# Ayyyyddd: Sampling date, year (yyyy), then day (ddd), between 1 and 365.
# h00v00: Identifies the grid tile (see https://lpdaac.usgs.gov/dataset_discovery/modis)
# vvv: Data version
# yyyydddhhmmss: Date data were processed, year, day, hour, minute, second.
dayOffset <- 121
dayRange <- 18:31
collectionDates <- paste("May", dayRange, "_2012", sep = "")
collectionDatesPOSIX <- as.POSIXct(paste("2012-05-", dayRange, sep = ""))
splitTemperaturesBySatellite <- lapply(c(Terra = "MOD11A1.A2012", Aqua = "MYD11A1.A2012"), function(searchString) {
temperatureFiles <- list.files(path = rawDataFilesLocation, pattern = searchString, full.names = TRUE)
subFiles <- sapply(paste("A2012", dayOffset + dayRange, sep = ""), grep, x = temperatureFiles, value = TRUE)
temperatures <- lapply(subFiles, MODIS::getSds)
splitTemperatures <- split(temperatures, f = factor(substr(subFiles, start = 0, stop = gregexpr(pattern = ".h2", text = subFiles[[1]])[[1]] - 1)))
names(splitTemperatures) <- collectionDates
splitTemperatures
})
indiaPolygons <- raster::getData(country = "IND", level = 2)
indiaPolygonsOtherCRS <- spTransform(indiaPolygons, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE)))
largeWesternMahaPolygonEdges <- rbind(c(21, 72.76), c(21, 76.5), c(17.0, 76.5), c(17.0, 72.6), c(21, 72.6))
largeWesternMahaPolygonEdges <- largeWesternMahaPolygonEdges[ , 2:1]
largeWesternMahaPolygon <- SpatialPolygons(Srl = list(Polygons(list(Polygon(coords = largeWesternMahaPolygonEdges)), ID = "Mumbai")))
crs(largeWesternMahaPolygon) <- crs(indiaPolygons)
largeWesternMahaPolygonOtherCRS <- spTransform(largeWesternMahaPolygon, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE)))
lastWeekIndices <- collectionDatesPOSIX >= as.POSIXct("2012-05-25")
indiaTemperaturesAndTimesLastWeek <- lapply(seq_along(splitTemperaturesBySatellite$Aqua)[lastWeekIndices], function(var1) {
aquaRasters <- funToCreateRaster(splitTemperaturesBySatellite$Aqua[[var1]], polygonBound = largeWesternMahaPolygonOtherCRS)
terraRasters <- funToCreateRaster(splitTemperaturesBySatellite$Terra[[var1]], polygonBound = largeWesternMahaPolygonOtherCRS)
if (sum(!is.na(values(aquaRasters$temperatureRaster))) >= sum(!is.na(values(terraRasters$temperatureRaster)))) {
cat("Returning Aqua!\n")
c(aquaRasters, satellite = "Aqua")
} else {
cat("Returning Terra!\n")
c(terraRasters, satellite = "Terra")
}
})
indiaTemperaturesLastWeek <- lapply(indiaTemperaturesAndTimesLastWeek, function(x) x$temperatureRaster)
names(indiaTemperaturesLastWeek) <- as.character(collectionDatesPOSIX[lastWeekIndices])
indiaTimesLastWeek <- lapply(indiaTemperaturesAndTimesLastWeek, function(x) x$hourRaster)
satellitePerDayLastWeek <- sapply(indiaTemperaturesAndTimesLastWeek, function(x) x$satellite)
citiesNames <- c("Mumbai City", "Pune", "Nashik", "Ahmadnagar", "Aurangabad")
cityPoints <- mapmisc::geocode(x = citiesNames)
cityPointsOtherCRS <- spTransform(cityPoints, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE)))
indiaExtent <- extent(c(xmin = min(largeWesternMahaPolygonEdges[ , 1]), xmax = max(largeWesternMahaPolygonEdges[ , 1]), ymin = min(largeWesternMahaPolygonEdges[ , 2]), ymax = max(largeWesternMahaPolygonEdges[ , 2])))
indiaRasterReprojected <- raster(x = indiaExtent, nrows = nrow(indiaTemperaturesLastWeek[[1]]), ncols = ncol(indiaTemperaturesLastWeek[[1]]), crs = crs(indiaPolygons))
# A bug appears when we use crs(cityPoints) instead of crs(indiaPolygons). The printout for both crs appears to be the same, but we get
# > identical(crs(indiaPolygons), crs(cityPoints))
# [1] FALSE
# This bug results probably from the update in gdal.
# We observe the following:
# > foo <- crs(cityPoints)
# > bar <- crs(indiaPolygons)
# > foo
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs
# > bar
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs
# > slotNames(foo)
# [1] "projargs"
# > slotNames(bar)
# [1] "projargs"
# > foo@projargs
# [1] "+proj=longlat +datum=WGS84 +no_defs"
# > bar@projargs
# [1] "+proj=longlat +datum=WGS84 +no_defs"
# > identical(bar@projargs, foo@projargs)
# [1] TRUE
# > object.size(foo)
# 1496 bytes
# > object.size(bar)
# 1560 bytes
# The only component of both objects are identical, but the objects have different size!
lapply(c("2012-05-27", "2012-05-28", "2012-05-29"), function(dateName) {
rasterReprojected <- raster::projectRaster(from = indiaTemperaturesLastWeek[[dateName]], to = indiaRasterReprojected)
ecol <- mapmisc::colourScale(values(rasterReprojected), col = "Spectral", breaks = c(20, 25, 30, 35, 40, 45, 50, 55, 60), style='fixed', rev = TRUE, opacity = 0.5)
filename <- paste("outputFiles/temperatures", dateName, ".jpg", sep = "")
jpeg(filename = filename, width = 1200, height = 1200)
raster::plot(rasterReprojected, legend = FALSE, cex.axis = 2.5, col = ecol$col, breaks = ecol$breaks)
plot(indiaPolygons, add = TRUE)
# raster::plot(indiaTemperaturesLastWeek[[dateName]], legend.only = TRUE, legend.width = 4, axis.args = list(cex.axis = 3))
mapmisc::legendBreaks("bottomleft", ecol, title = NULL, bg = "white", cex = 2.5, bty = "n")
plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE)
text(x = cityPoints@coords[ , 1], y = cityPoints@coords[ , 2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 4, pos = 4)
scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3)
dev.off()
NULL
})
We now produce the plots used to illustrate the validation analysis following the application we described. Uncommenting the lines used in the creation of the plots will add city names to the graphs.
smallWesternMahaPolygonEdges <- rbind(c(19.55, 72.76), c(19.55, 74), c(18.35, 74), c(18.35, 72.76), c(19.55, 72.76))
smallWesternMahaPolygonEdges <- smallWesternMahaPolygonEdges[ , 2:1]
smallWesternMahaPolygon <- SpatialPolygons(Srl = list(Polygons(list(Polygon(coords = smallWesternMahaPolygonEdges)), ID = "Mumbai")))
crs(smallWesternMahaPolygon) <- crs(indiaPolygons)
smallWesternMahaPolygonOtherCRS <- spTransform(smallWesternMahaPolygon, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE)))
secondToLastWeekIndices <- (collectionDatesPOSIX >= as.POSIXct("2012-05-18")) & (collectionDatesPOSIX < as.POSIXct("2012-05-25"))
indiaTemperaturesAndTimesSecondToLastWeek <- lapply(seq_along(splitTemperaturesBySatellite$Aqua)[secondToLastWeekIndices], function(var1) {
aquaRasters <- funToCreateRaster(splitTemperaturesBySatellite$Aqua[[var1]], polygonBound = smallWesternMahaPolygonOtherCRS)
terraRasters <- funToCreateRaster(splitTemperaturesBySatellite$Terra[[var1]], polygonBound = smallWesternMahaPolygonOtherCRS)
if (sum(!is.na(values(aquaRasters$temperatureRaster))) >= sum(!is.na(values(terraRasters$temperatureRaster)))) {
cat("Returning Aqua!\n")
c(aquaRasters, satellite = "Aqua")
} else {
cat("Returning Terra!\n")
c(terraRasters, satellite = "Terra")
}
})
names(indiaTemperaturesAndTimesSecondToLastWeek) <- collectionDatesPOSIX[secondToLastWeekIndices]
indiaTemperaturesSecondToLastWeek <- lapply(indiaTemperaturesAndTimesSecondToLastWeek, function(x) x$temperatureRaster)
names(indiaTemperaturesSecondToLastWeek) <- names(indiaTemperaturesAndTimesSecondToLastWeek)
indiaTimesSecondToLastWeek <- lapply(indiaTemperaturesAndTimesSecondToLastWeek, function(x) x$hourRaster)
satellitePerDaySecondToLastWeek <- sapply(indiaTemperaturesAndTimesSecondToLastWeek, function(x) x$satellite)
smallRasterExtent <- extent(c(xmin = min(smallWesternMahaPolygonEdges[ , 1]), xmax = max(smallWesternMahaPolygonEdges[ , 1]), ymin = min(smallWesternMahaPolygonEdges[ , 2]), ymax = max(smallWesternMahaPolygonEdges[ , 2])))
emptyRaster <- raster(x = smallRasterExtent, nrows = 120, ncols = 120, crs = crs(indiaPolygons))
######## Producing plot to represent validation test set ########
may21raster <- may21rasterAddedMissing <- indiaTemperaturesSecondToLastWeek[["2012-05-21"]]
may28raster <- crop(x = indiaTemperaturesLastWeek[["2012-05-28"]], y = extent(may21raster))
indicesForKnownRemovedValues <- which(is.na(values(may28raster)) & !is.na(values(may21rasterAddedMissing)))
values(may21rasterAddedMissing) <- replace(values(may21rasterAddedMissing), indicesForKnownRemovedValues, NA)
may21rasterAddedMissingReproj <- projectRaster(from = may21rasterAddedMissing, to = emptyRaster)
may21rasterReproj <- projectRaster(from = may21raster, to = emptyRaster)
ecolBasicValid <- mapmisc::colourScale(values(may21rasterReproj), col = "Spectral", dec = 0, breaks = c(20, 25, 30, 35, 40, 45, 50, 55, 60), style = 'fixed', rev = TRUE, opacity = 0.5)
jpeg("outputFiles/may21rasterOriginal.jpg", width = 1200, height = 1200)
raster::plot(may21rasterReproj, legend = FALSE, cex.axis = 2.5, col = ecolBasicValid$col, breaks = ecolBasicValid$breaks)
raster::plot(indiaPolygons, add = TRUE)
# raster::plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE)
# text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 4, pos = 4)
mapmisc::legendBreaks("bottomleft", ecolBasicValid, title = NULL, bg = "white", cex = 2.5, bty = "n")
scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3)
dev.off()
jpeg("outputFiles/may21rasterAddedMissing.jpg", width = 1200, height = 1200)
raster::plot(may21rasterAddedMissingReproj, legend = FALSE, cex.axis = 2.5, col = ecolBasicValid$col, breaks = ecolBasicValid$breaks)
raster::plot(indiaPolygons, add = TRUE)
# raster::plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE)
# text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 4, pos = 4)
mapmisc::legendBreaks("bottomleft", ecolBasicValid, title = NULL, bg = "white", cex = 2.5, bty = "n")
scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3)
dev.off()
Several prediction error summaries can be found in the text. The next chunk contains the code used to calculate them.
temperaturesMay21 <- indiaTemperaturesSecondToLastWeek[["2012-05-21"]]
temperaturesMay21reproj <- raster(x = extent(testDataValidation@sp), nrows = nrow(temperaturesMay21), ncols = ncol(temperaturesMay21), crs = crs(indiaPolygonsOtherCRS))
temperaturesMay21reproj <- raster::projectRaster(from = temperaturesMay21, to = temperaturesMay21reproj)
recordedTemperaturesInMissingZone <- raster::extract(x = temperaturesMay21reproj, y = testDataValidation@sp@coords)
differences <- indiaAnalysisValidation$predMoments$Mean - recordedTemperaturesInMissingZone
MSPE <- mean(differences^2, na.rm = TRUE)
MedSPE <- median(differences^2, na.rm = TRUE)
quantile(abs(differences), probs = 0.9)
# What is the coordinate of the tile that produced the largest absolute difference?
testDataSpReproj <- spTransform(testDataValidation@sp, CRSobj = crs(indiaPolygons))
testDataSpReproj@coords[which.max(abs(differences)),]
spObjectReprojected <- sp::spTransform(testDataValidation@sp, CRSobj = crs(indiaPolygons))
spObjectReprojected@coords[which.min(indiaAnalysisValidation$predMoments$Mean), ]
fieldValues <- differences[!is.na(differences)]
sqDiffRaster <- raster::rasterize(x = spObjectReprojected, y = emptyRaster, field = fieldValues^2)
diffRaster <- raster::rasterize(x = spObjectReprojected, y = emptyRaster, field = fieldValues)
The previous results were all for IS-MRA. The next chunk evaluates the predictive performance of INLA-SPDE. We re-create the mesh and transformation matrix used when we fitted the model.
library(INLA)
trainingDataValidation@sp <- sp::spTransform(trainingDataValidation@sp, CRSobj = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
testDataValidation@sp <- sp::spTransform(testDataValidation@sp, CRSobj = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
timeVecTraining <- (as.numeric(time(trainingDataValidation)) - min(as.numeric(time(trainingDataValidation))))/(3600*24) + 1
timeVecTest <- (as.numeric(time(testDataValidation)) - min(as.numeric(time(trainingDataValidation))))/(3600*24) + 1
knots <- seq(1, max(timeVecTraining), length = max(timeVecTraining))
mesh1 <- inla.mesh.1d(loc = knots, degree = 2, boundary = "free")
## generate space mesh
mesh2 <- inla.mesh.2d(loc = trainingDataValidation@sp@coords[timeVecTraining == 1, ], cutoff = 0.01, offset = c(0.1, 0.2), max.n = 2000)
# range0 and sigma0 control the prior means for the range and scale parameters.
# See Lindgren INLA tutorial page 5.
d <- 1
alpha <- 2
kappa <- 1 # = 1/(range parameter in my model)
spatialSmoothness <- alpha - d/2 # cf p.3 INLA tutorial
loghyperparaSDinMyModel <- log(10)
# range0 and sigma0 seem to be the prior means...
range0 <- sqrt(8 * spatialSmoothness)/kappa # sqrt(8 * spatial smoothness) / Kappa. In my model, I use 1 as prior mean for spatial range and fix smoothness at 1.5. This means Kappa = 1.
sigma0 <- 1
lkappa0 <- log(8 * spatialSmoothness)/2 - log(range0)
ltau0 <- 0.5*log(gamma(spatialSmoothness)/(gamma(alpha)*(4*pi)^(d/2))) - log(sigma0) - spatialSmoothness * lkappa0
## build the spatial spde
spde <- inla.spde2.matern(mesh2, B.tau = matrix(c(ltau0, -1, spatialSmoothness), 1, 3),
B.kappa = matrix(c(lkappa0, 0, -1), 1,3),
theta.prior.mean = c(0,0), theta.prior.prec = c(1/loghyperparaSDinMyModel^2, 1/loghyperparaSDinMyModel^2))
## build the space time indices
STindex <- inla.spde.make.index("space", n.spde = spde$n.spde, n.group = mesh1$m)
## Link data and process
Atraining <- inla.spde.make.A(mesh2, loc = trainingDataValidation@sp@coords, group = timeVecTraining, group.mesh = mesh1)
Atest <- inla.spde.make.A(mesh2, loc = testDataValidation@sp@coords, group = timeVecTest, group.mesh = mesh1)
stackTraining <- inla.stack(data = list(y = trainingDataValidation@data$y), A = list(Atraining, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
effects = list(
c(STindex, list(intercept = 1)),
list(landCover3 = trainingDataValidation@data$landCover2), # The mismatch between the element name and the name of the extracted column is to reflect an error in the original fitting, where landCover2 was encoded as landCover3.
list(landCover4 = trainingDataValidation@data$landCover4),
list(landCover5 = trainingDataValidation@data$landCover5),
list(landCover8 = trainingDataValidation@data$landCover8),
list(landCover9 = trainingDataValidation@data$landCover9),
list(landCover10 = trainingDataValidation@data$landCover10),
list(landCover11 = trainingDataValidation@data$landCover11),
list(landCover12 = trainingDataValidation@data$landCover12),
list(landCover13 = trainingDataValidation@data$landCover13),
list(landCover14 = trainingDataValidation@data$landCover14),
list(landCover15 = trainingDataValidation@data$landCover15),
list(elevation = trainingDataValidation@data$elevation),
list(Aqua = trainingDataValidation@data$Aqua),
list(time2 = trainingDataValidation@data$time2),
list(time3 = trainingDataValidation@data$time3),
list(time4 = trainingDataValidation@data$time4),
list(time5 = trainingDataValidation@data$time5),
list(time6 = trainingDataValidation@data$time6),
list(time7 = trainingDataValidation@data$time7)
), tag = "est")
stackTest <- inla.stack(
data = list(y = NA),
A = list(Atest, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 , 1, 1, 1, 1),
effects = list(
c(STindex, list(intercept = 1)),
list(landCover3 = testDataValidation@data$landCover2), # The mismatch between the element name and the name of the extracted column is to reflect an error in the original fitting, where landCover2 was encoded as landCover3.
list(landCover4 = testDataValidation@data$landCover4),
list(landCover5 = testDataValidation@data$landCover5),
list(landCover8 = testDataValidation@data$landCover8),
list(landCover9 = testDataValidation@data$landCover9),
list(landCover10 = testDataValidation@data$landCover10),
list(landCover11 = testDataValidation@data$landCover11),
list(landCover12 = testDataValidation@data$landCover12),
list(landCover13 = testDataValidation@data$landCover13),
list(landCover14 = testDataValidation@data$landCover14),
list(landCover15 = testDataValidation@data$landCover15),
list(elevation = testDataValidation@data$elevation),
list(Aqua = testDataValidation@data$Aqua),
list(time2 = testDataValidation@data$time2),
list(time3 = testDataValidation@data$time3),
list(time4 = testDataValidation@data$time4),
list(time5 = testDataValidation@data$time5),
list(time6 = testDataValidation@data$time6),
list(time7 = testDataValidation@data$time7)
),
tag = 'predictions')
combinedStack <- inla.stack(stackTraining, stackTest)
stackIndex <- inla.stack.index(combinedStack, "predictions")$data
preds <- SPDEresult$summary.linear.predictor
differencesSPDE <- preds$mean[stackIndex] - recordedTemperaturesInMissingZone
MSPE_SPDE <- mean(differencesSPDE^2, na.rm = TRUE)
MedSPE_SPDE <- median(differencesSPDE^2, na.rm = TRUE)
smallRasterExtent <- raster::extent(c(xmin = min(smallWesternMahaPolygonEdges[ , 1]), xmax = max(smallWesternMahaPolygonEdges[ , 1]), ymin = min(smallWesternMahaPolygonEdges[ , 2]), ymax = max(smallWesternMahaPolygonEdges[ , 2])))
fieldValuesSPDE <- differencesSPDE[!is.na(differencesSPDE)]
spObjectReprojectedSPDE <- sp::spTransform(testDataValidation@sp[!is.na(differencesSPDE)], CRSobj = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
emptyRaster <- raster(x = smallRasterExtent, nrows = 120, ncols = 120, crs = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
diffRasterSPDE <- raster::rasterize(x = spObjectReprojectedSPDE, y = emptyRaster, field = fieldValuesSPDE)
ecolSPDE <- mapmisc::colourScale(values(diffRasterSPDE), col = "Spectral", dec = 0, breaks = 10, rev = TRUE, style = "equal", opacity = 0.5)
ecolSPDE <- ecolSPDE[c("col", "breaks")]
ecolINLAMRA <- mapmisc::colourScale(values(diffRaster), col = "Spectral", dec = 0, breaks = 10, rev = TRUE, style = "equal", opacity = 0.5)
ecolINLAMRA <- ecolINLAMRA[c("col", "breaks")]
newBoundaries <- range(c(ecolINLAMRA$breaks, ecolSPDE$breaks))
combinedBreaks <- seq(from = newBoundaries[[1]], to = newBoundaries[[2]], length.out = length(ecolINLAMRA$breaks))
jpeg(file = "outputFiles/predErrorOnMay21withTimeCovarSPDE.jpg", width = 1200, height = 1200)
raster::plot(diffRasterSPDE, legend = FALSE, cex.axis = 2.5, col = ecolINLAMRA$col, breaks = combinedBreaks)
plot(indiaPolygons, add = TRUE)
# plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE)
# text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = citiesNames, offset = 2, cex = 5, pos = 4)
mapmisc::legendBreaks("bottomleft", list(breaks = round(combinedBreaks, digits = 1), col = ecolINLAMRA$col), title = NULL, bg = "white", cex = 2.5, bty = "n")
mapmisc::scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3)
dev.off()
jpeg(file = "outputFiles/predErrorOnMay21withTimeCovar.jpg", width = 1200, height = 1200)
raster::plot(diffRaster, legend = FALSE, cex.axis = 2.5, col = ecolINLAMRA$col, breaks = combinedBreaks)
plot(indiaPolygons, add = TRUE)
mapmisc::legendBreaks("bottomleft", list(breaks = round(combinedBreaks, digits = 1), col = ecolINLAMRA$col), title = NULL, bg = "white", cex = 2.5, bty = "n")
# raster::plot(diffRaster, legend.only = TRUE, legend.width = 4, axis.args = list(cex.axis = 3))
# plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE)
# text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 5, pos = 4)
mapmisc::scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3)
dev.off()
The last graph we produced presents prediction errors for IS-MRA across the entire prediction region. We could not produce it before because we wanted to use the same colour scale as for the analogous graph for SPDE, to make them easily comparable when read side-by-side. We now move on to the main analysis, on the large dataset.
spAndTimeCorr <- lapply(list(space = list(paraName = "space.rho", smoothness = 1.5, distances = c(close = 1, far = 10)), time = list(paraName = "time.rho", smoothness = 0.5, distances = c(close = 1, far = 7))), function(component) {
sapply(component$distances, function(distanceMeasure) {
MRAinla::maternCov(d = distanceMeasure, rho = exp(indiaAnalysisMain$hyperMarginalMoments[component$paraName, "Mean"]), smoothness = component$smoothness, scale = 1)
})
})
######### Preparing parameter and hyperparameter moments table ########
hyperMoments <- subset(indiaAnalysisMain$hyperMarginalMoments[c("space.rho", "time.rho", "scale"), ], select = -Skewness)
## Delta method to restore the original scale
hyperMoments$StdDev <- hyperMoments$StdDev * exp(hyperMoments$Mean)
# mean(rho) = exp(mean(log rho) + 0.5 var(log rho))
hyperMoments$Mean <- exp(hyperMoments$Mean + 0.5 * hyperMoments$StdDev^2)
hyperMoments[ , "CredInt_2.5%"] <- exp(hyperMoments[ , "CredInt_2.5%"])
hyperMoments[ , "CredInt_97.5%"] <- exp(hyperMoments[ , "CredInt_97.5%"])
covariateMoments <- indiaAnalysisMain$FEmarginalMoments
combinedMoments <- rbind(covariateMoments, hyperMoments)
landCoverNames <- c("Water", "Evergreen Needleleaf Forests", "Evergreen Broadleaf Forests", "Deciduous Needleleaf Forests", "Deciduous Broadleaf Forests", "Mixed Forests", "Closed Shrublands", "Open Shrublands", "Woody Savannas", "Savannas", "Grasslands", "Permanent Wetlands", "Croplands", "Urban and Built-up Lands", "Cropland/Natural Vegetation Mosaics", "Non-Vegetated Lands", "Unclassified")
names(landCoverNames) <- paste("landCover", c(0:15, 255), sep = "")
landCoverPositions <- grep(pattern = "landCover", x = rownames(combinedMoments))
landCoverValues <- grep(pattern = "landCover", x = rownames(combinedMoments), value = TRUE)
rownames(combinedMoments) <- replace(
rownames(combinedMoments),
landCoverPositions,
unname(landCoverNames[landCoverValues])
)
rownames(combinedMoments) <- replace(
rownames(combinedMoments),
match("elevation", rownames(combinedMoments)),
"Elevation"
)
rownames(combinedMoments) <- replace(
rownames(combinedMoments),
match("Aqua", rownames(combinedMoments)),
"Satellite: Aqua"
)
datesCovered <- paste("May", 26:31)
names(datesCovered) <- paste("time", 2:7, sep = "")
rownames(combinedMoments) <- replace(
rownames(combinedMoments),
grep(pattern = "time[0-9]$", x = rownames(combinedMoments)),
unname(datesCovered[grep(pattern = "time[0-9]$", x = rownames(combinedMoments), value = TRUE)])
)
reorderedCombined <- combinedMoments[c("Intercept", "Elevation", "Satellite: Aqua", intersect(rownames(combinedMoments), landCoverNames), intersect(rownames(combinedMoments), datesCovered), rownames(hyperMoments)), ]
latexCode <- Hmisc::latex(Hmisc::format.df(reorderedCombined, dec = 3), rgroup = c("", "Land cover", "Time", "Hyperparameters"), n.rgroup = c(3, 14, 6, 3), file = "outputFiles/JASAtable1.tex")
######## END: Preparing parameter and hyperparameter moments table ########
table((indiaAnalysisMain$predMoments$Mean - min(trainingDataMain@data[, "y"])) < -1)
extremeValueIndices <- which(indiaAnalysisMain$predMoments$Mean < (min(trainingDataMain@data[, "y"]) - 1))
mostExtremeIndex <- which.min(indiaAnalysisMain$predMoments$Mean)
predictionDataMain@data[mostExtremeIndex, ]
predictionOrder <- order(indiaAnalysisMain$predMoments$Mean)
predictionDataMain@sp@coords[predictionOrder[1:2],]
extremeDataset <- predictionDataMain[extremeValueIndices]
extremeDataset@data <- data.frame(Temperature = indiaAnalysisMain$predMoments$Mean[extremeValueIndices])
# jpeg("outputFiles/extremeValuePositions.jpeg", width = 1000, height = 1000)
# stplot(extremeDataset)
# dev.off()
######## Plotting predictions ########
indiaPolygons <- raster::getData(country = "IND", level = 2)
trainingDataMainReproject <- trainingDataMain
trainingDataMainReproject@sp <- sp::spTransform(trainingDataMain@sp, crs(indiaPolygons))
predictionDataMainReproject <- predictionDataMain
predictionDataMainReproject@sp <- sp::spTransform(predictionDataMain@sp, crs(indiaPolygons))
# Predicted values outside the range of observed values will be omitted to get a more appealing visualisation. The following lines compute the number of omitted values:
table(indiaAnalysisMain$predMoments$Mean[which(indiaAnalysisMain$predData$time == as.POSIXct("2012-05-28"))] <= min(indiaAnalysisMain$data$spObject@data$y[which(indiaAnalysisMain$data$time == as.POSIXct("2012-05-28"))]))
table(indiaAnalysisMain$predMoments$Mean[which(indiaAnalysisMain$predData$time == as.POSIXct("2012-05-28"))] >= max(indiaAnalysisMain$data$spObject@data$y[which(indiaAnalysisMain$data$time == as.POSIXct("2012-05-28"))]))
plot.INLAMRA(indiaAnalysisMain,
filename = "outputFiles/predictionsIndiaDataMay28mainAnalysisJoint_Larger.jpg",
type = "joint",
polygonsToOverlay = indiaPolygons,
control = list(
graphicsEngine = jpeg,
controlForScaleBar = list(
pos = "topleft",
cex = 2,
pt.cex = 1.5,
seg.len = 3
),
controlForRasterLegend = list(
pos = "bottomleft",
title = NULL,
bg = "white",
cex = 2.5
),
controlForRasterColourScale = list(
col = "Spectral",
breaks = seq(from = 20, to = 60, by = 5), style='fixed',
rev = TRUE,
dec = 1,
opacity = 0.5
),
controlForRasterPlot = list(
cex.axis = 2,
cex.main = 3,
main = NULL
),
resolutionInMeters = 1000,
trim = 2,
timesToPlot = unique(time(predictionDataMain))
),
width = 1100,
height = 1200
)
# There was a bug in the estimation of the SDs in an early version of the software: to get the real value, we must apply the following transformation:
# sqrt(indiaAnalysisMain$predMoments$SD - errorSD^2) + errorSD
errorSD <- 0.5 # cf. vignette "ISMRAmainAnalysis_Script.Rmd"
indiaAnalysisMainCopy <- indiaAnalysisMain
indiaAnalysisMainCopy$predMoments$SD <- sqrt(indiaAnalysisMain$predMoments$SD - errorSD^2) + errorSD
plot.INLAMRA(indiaAnalysisMainCopy,
filename = "outputFiles/predictionsIndiaDataMay28mainAnalysisSDs_LargerFixed.jpg",
type = "SD",
polygonsToOverlay = indiaPolygons,
control = list(
graphicsEngine = jpeg,
controlForScaleBar = list(
pos = "topleft",
cex = 2,
pt.cex = 1.5,
seg.len = 3
),
controlForRasterLegend = list(
pos = "bottomleft",
title = NULL,
bg = "white",
cex = 2.5
),
controlForRasterColourScale = list(
col = "Spectral",
breaks = c(0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5), style='fixed',
rev = TRUE,
dec = 1,
opacity = 1
),
controlForRasterPlot = list(
cex.axis = 2,
cex.main = 3,
main = NULL
),
resolutionInMeters = 1000,
timesToPlot = unique(time(predictionDataMain))
),
width = 1100,
height = 1200
)
plot.INLAMRA(indiaAnalysisMain,
filename = "outputFiles/predictionsIndiaDataMay28mainAnalysisTraining_Larger.jpg",
type = "training",
polygonsToOverlay = indiaPolygons,
control = list(
graphicsEngine = jpeg,
controlForScaleBar = list(
pos = "topleft",
cex = 2,
pt.cex = 1.5,
seg.len = 3
),
controlForRasterLegend = list(
pos = "bottomleft",
title = NULL,
bg = "white",
cex = 2.5
),
controlForRasterColourScale = list(
col = "Spectral",
breaks = seq(from = 20, to = 60, by = 5), style='fixed',
rev = TRUE,
dec = 0,
opacity = 0.5
),
controlForRasterPlot = list(
cex.axis = 2,
cex.main = 3,
main = NULL
),
resolutionInMeters = 1000,
timesToPlot = unique(time(predictionDataMain))
),
width = 1100,
height = 1200
)
villandre/ISMRAanalyses documentation built on Dec. 23, 2021, 3:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
We start by loading the required libraries.
if ("package:ISMRAanalyses" %in% search()) { detach("package:ISMRAanalyses", unload = TRUE) } library(ISMRAanalyses) library(sp) library(raster) library(MODIS) library(rgdal) library(rgeos) library(spacetime) library(MRAinla) library(RhpcBLASctl) library(geoR) library(maptools) library(doParallel) library(mapmisc) library(xtable) # Input the correct working directory for your system # **The working directory should have two subfolders: "data" for datasets and "outputFiles" for model fitting outputs and graphs.** setwd("/home/luc/INLAMRAfiles/INLAMRApaper1/realData") RandomFields::RFoptions(cores = 1) blas_set_num_threads(1) omp_set_num_threads(1)
Now, we load all datasets required for producing the results presented in the article.
# Input in variable rawDataFilesLocation the folder name where MODIS data files are found rawDataFilesLocation <- "/store/luc/rawDataFiles" trainingDataValidationName <- load("data/mainDataCompleteMap_May18_24.Rdata") trainingDataValidation <- get(trainingDataValidationName) rm(list = trainingDataValidationName) testDataValidationName <- load("data/testDataMay21_May18_24.Rdata") testDataValidation <- get(testDataValidationName) rm(list = testDataValidationName) indiaAnalysisValidationName <- load("outputFiles/INLAMRA_validationAnalysis.Rdata") indiaAnalysisValidation <- get(indiaAnalysisValidationName) if (indiaAnalysisValidationName != "indiaAnalysisValidation") { rm(list = indiaAnalysisValidationName) } SPDEresultName <- load("outputFiles/inlaFitForValidationAnalysis.Rdata") SPDEresult <- get(SPDEresultName) rm(list = SPDEresultName) predictionDataMainName <- load("data/testDataMay28_May25_31_Larger.Rdata") predictionDataMain <- get(predictionDataMainName) rm(list = predictionDataMainName) trainingDataMainName <- load("data/mainDataCompleteMap_May25_31_Larger.Rdata") trainingDataMain <- get(trainingDataMainName) rm(list = trainingDataMainName) indiaAnalysisMainName <- load("outputFiles/indiaAnalysisMay28evenLarger_May25_31_VersionApril13.Rdata") indiaAnalysisMain <- get(indiaAnalysisMainName) rm(list = indiaAnalysisMainName)
We now produce the basic raster plots for the observed land surface temperatures, with cities and geopolitical boundaries overlaid.
# Naming convention: nnnnnnn.Ayyyyddd.h00v00.vvv.yyyydddhhmmss. # nnnnnnn: Product name # Ayyyyddd: Sampling date, year (yyyy), then day (ddd), between 1 and 365. # h00v00: Identifies the grid tile (see https://lpdaac.usgs.gov/dataset_discovery/modis) # vvv: Data version # yyyydddhhmmss: Date data were processed, year, day, hour, minute, second. dayOffset <- 121 dayRange <- 18:31 collectionDates <- paste("May", dayRange, "_2012", sep = "") collectionDatesPOSIX <- as.POSIXct(paste("2012-05-", dayRange, sep = "")) splitTemperaturesBySatellite <- lapply(c(Terra = "MOD11A1.A2012", Aqua = "MYD11A1.A2012"), function(searchString) { temperatureFiles <- list.files(path = rawDataFilesLocation, pattern = searchString, full.names = TRUE) subFiles <- sapply(paste("A2012", dayOffset + dayRange, sep = ""), grep, x = temperatureFiles, value = TRUE) temperatures <- lapply(subFiles, MODIS::getSds) splitTemperatures <- split(temperatures, f = factor(substr(subFiles, start = 0, stop = gregexpr(pattern = ".h2", text = subFiles[[1]])[[1]] - 1))) names(splitTemperatures) <- collectionDates splitTemperatures }) indiaPolygons <- raster::getData(country = "IND", level = 2) indiaPolygonsOtherCRS <- spTransform(indiaPolygons, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE))) largeWesternMahaPolygonEdges <- rbind(c(21, 72.76), c(21, 76.5), c(17.0, 76.5), c(17.0, 72.6), c(21, 72.6)) largeWesternMahaPolygonEdges <- largeWesternMahaPolygonEdges[ , 2:1] largeWesternMahaPolygon <- SpatialPolygons(Srl = list(Polygons(list(Polygon(coords = largeWesternMahaPolygonEdges)), ID = "Mumbai"))) crs(largeWesternMahaPolygon) <- crs(indiaPolygons) largeWesternMahaPolygonOtherCRS <- spTransform(largeWesternMahaPolygon, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE))) lastWeekIndices <- collectionDatesPOSIX >= as.POSIXct("2012-05-25") indiaTemperaturesAndTimesLastWeek <- lapply(seq_along(splitTemperaturesBySatellite$Aqua)[lastWeekIndices], function(var1) { aquaRasters <- funToCreateRaster(splitTemperaturesBySatellite$Aqua[[var1]], polygonBound = largeWesternMahaPolygonOtherCRS) terraRasters <- funToCreateRaster(splitTemperaturesBySatellite$Terra[[var1]], polygonBound = largeWesternMahaPolygonOtherCRS) if (sum(!is.na(values(aquaRasters$temperatureRaster))) >= sum(!is.na(values(terraRasters$temperatureRaster)))) { cat("Returning Aqua!\n") c(aquaRasters, satellite = "Aqua") } else { cat("Returning Terra!\n") c(terraRasters, satellite = "Terra") } }) indiaTemperaturesLastWeek <- lapply(indiaTemperaturesAndTimesLastWeek, function(x) x$temperatureRaster) names(indiaTemperaturesLastWeek) <- as.character(collectionDatesPOSIX[lastWeekIndices]) indiaTimesLastWeek <- lapply(indiaTemperaturesAndTimesLastWeek, function(x) x$hourRaster) satellitePerDayLastWeek <- sapply(indiaTemperaturesAndTimesLastWeek, function(x) x$satellite) citiesNames <- c("Mumbai City", "Pune", "Nashik", "Ahmadnagar", "Aurangabad") cityPoints <- mapmisc::geocode(x = citiesNames) cityPointsOtherCRS <- spTransform(cityPoints, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE))) indiaExtent <- extent(c(xmin = min(largeWesternMahaPolygonEdges[ , 1]), xmax = max(largeWesternMahaPolygonEdges[ , 1]), ymin = min(largeWesternMahaPolygonEdges[ , 2]), ymax = max(largeWesternMahaPolygonEdges[ , 2]))) indiaRasterReprojected <- raster(x = indiaExtent, nrows = nrow(indiaTemperaturesLastWeek[[1]]), ncols = ncol(indiaTemperaturesLastWeek[[1]]), crs = crs(indiaPolygons)) # A bug appears when we use crs(cityPoints) instead of crs(indiaPolygons). The printout for both crs appears to be the same, but we get # > identical(crs(indiaPolygons), crs(cityPoints)) # [1] FALSE # This bug results probably from the update in gdal. # We observe the following: # > foo <- crs(cityPoints) # > bar <- crs(indiaPolygons) # > foo # CRS arguments: +proj=longlat +datum=WGS84 +no_defs # > bar # CRS arguments: +proj=longlat +datum=WGS84 +no_defs # > slotNames(foo) # [1] "projargs" # > slotNames(bar) # [1] "projargs" # > foo@projargs # [1] "+proj=longlat +datum=WGS84 +no_defs" # > bar@projargs # [1] "+proj=longlat +datum=WGS84 +no_defs" # > identical(bar@projargs, foo@projargs) # [1] TRUE # > object.size(foo) # 1496 bytes # > object.size(bar) # 1560 bytes # The only component of both objects are identical, but the objects have different size! lapply(c("2012-05-27", "2012-05-28", "2012-05-29"), function(dateName) { rasterReprojected <- raster::projectRaster(from = indiaTemperaturesLastWeek[[dateName]], to = indiaRasterReprojected) ecol <- mapmisc::colourScale(values(rasterReprojected), col = "Spectral", breaks = c(20, 25, 30, 35, 40, 45, 50, 55, 60), style='fixed', rev = TRUE, opacity = 0.5) filename <- paste("outputFiles/temperatures", dateName, ".jpg", sep = "") jpeg(filename = filename, width = 1200, height = 1200) raster::plot(rasterReprojected, legend = FALSE, cex.axis = 2.5, col = ecol$col, breaks = ecol$breaks) plot(indiaPolygons, add = TRUE) # raster::plot(indiaTemperaturesLastWeek[[dateName]], legend.only = TRUE, legend.width = 4, axis.args = list(cex.axis = 3)) mapmisc::legendBreaks("bottomleft", ecol, title = NULL, bg = "white", cex = 2.5, bty = "n") plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE) text(x = cityPoints@coords[ , 1], y = cityPoints@coords[ , 2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 4, pos = 4) scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3) dev.off() NULL })
We now produce the plots used to illustrate the validation analysis following the application we described. Uncommenting the lines used in the creation of the plots will add city names to the graphs.
smallWesternMahaPolygonEdges <- rbind(c(19.55, 72.76), c(19.55, 74), c(18.35, 74), c(18.35, 72.76), c(19.55, 72.76)) smallWesternMahaPolygonEdges <- smallWesternMahaPolygonEdges[ , 2:1] smallWesternMahaPolygon <- SpatialPolygons(Srl = list(Polygons(list(Polygon(coords = smallWesternMahaPolygonEdges)), ID = "Mumbai"))) crs(smallWesternMahaPolygon) <- crs(indiaPolygons) smallWesternMahaPolygonOtherCRS <- spTransform(smallWesternMahaPolygon, CRSobj = crs(readGDAL(splitTemperaturesBySatellite$Aqua[[1]][[1]]$SDS4gdal[1], as.is = TRUE))) secondToLastWeekIndices <- (collectionDatesPOSIX >= as.POSIXct("2012-05-18")) & (collectionDatesPOSIX < as.POSIXct("2012-05-25")) indiaTemperaturesAndTimesSecondToLastWeek <- lapply(seq_along(splitTemperaturesBySatellite$Aqua)[secondToLastWeekIndices], function(var1) { aquaRasters <- funToCreateRaster(splitTemperaturesBySatellite$Aqua[[var1]], polygonBound = smallWesternMahaPolygonOtherCRS) terraRasters <- funToCreateRaster(splitTemperaturesBySatellite$Terra[[var1]], polygonBound = smallWesternMahaPolygonOtherCRS) if (sum(!is.na(values(aquaRasters$temperatureRaster))) >= sum(!is.na(values(terraRasters$temperatureRaster)))) { cat("Returning Aqua!\n") c(aquaRasters, satellite = "Aqua") } else { cat("Returning Terra!\n") c(terraRasters, satellite = "Terra") } }) names(indiaTemperaturesAndTimesSecondToLastWeek) <- collectionDatesPOSIX[secondToLastWeekIndices] indiaTemperaturesSecondToLastWeek <- lapply(indiaTemperaturesAndTimesSecondToLastWeek, function(x) x$temperatureRaster) names(indiaTemperaturesSecondToLastWeek) <- names(indiaTemperaturesAndTimesSecondToLastWeek) indiaTimesSecondToLastWeek <- lapply(indiaTemperaturesAndTimesSecondToLastWeek, function(x) x$hourRaster) satellitePerDaySecondToLastWeek <- sapply(indiaTemperaturesAndTimesSecondToLastWeek, function(x) x$satellite) smallRasterExtent <- extent(c(xmin = min(smallWesternMahaPolygonEdges[ , 1]), xmax = max(smallWesternMahaPolygonEdges[ , 1]), ymin = min(smallWesternMahaPolygonEdges[ , 2]), ymax = max(smallWesternMahaPolygonEdges[ , 2]))) emptyRaster <- raster(x = smallRasterExtent, nrows = 120, ncols = 120, crs = crs(indiaPolygons)) ######## Producing plot to represent validation test set ######## may21raster <- may21rasterAddedMissing <- indiaTemperaturesSecondToLastWeek[["2012-05-21"]] may28raster <- crop(x = indiaTemperaturesLastWeek[["2012-05-28"]], y = extent(may21raster)) indicesForKnownRemovedValues <- which(is.na(values(may28raster)) & !is.na(values(may21rasterAddedMissing))) values(may21rasterAddedMissing) <- replace(values(may21rasterAddedMissing), indicesForKnownRemovedValues, NA) may21rasterAddedMissingReproj <- projectRaster(from = may21rasterAddedMissing, to = emptyRaster) may21rasterReproj <- projectRaster(from = may21raster, to = emptyRaster) ecolBasicValid <- mapmisc::colourScale(values(may21rasterReproj), col = "Spectral", dec = 0, breaks = c(20, 25, 30, 35, 40, 45, 50, 55, 60), style = 'fixed', rev = TRUE, opacity = 0.5) jpeg("outputFiles/may21rasterOriginal.jpg", width = 1200, height = 1200) raster::plot(may21rasterReproj, legend = FALSE, cex.axis = 2.5, col = ecolBasicValid$col, breaks = ecolBasicValid$breaks) raster::plot(indiaPolygons, add = TRUE) # raster::plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE) # text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 4, pos = 4) mapmisc::legendBreaks("bottomleft", ecolBasicValid, title = NULL, bg = "white", cex = 2.5, bty = "n") scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3) dev.off() jpeg("outputFiles/may21rasterAddedMissing.jpg", width = 1200, height = 1200) raster::plot(may21rasterAddedMissingReproj, legend = FALSE, cex.axis = 2.5, col = ecolBasicValid$col, breaks = ecolBasicValid$breaks) raster::plot(indiaPolygons, add = TRUE) # raster::plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE) # text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 4, pos = 4) mapmisc::legendBreaks("bottomleft", ecolBasicValid, title = NULL, bg = "white", cex = 2.5, bty = "n") scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3) dev.off()
Several prediction error summaries can be found in the text. The next chunk contains the code used to calculate them.
temperaturesMay21 <- indiaTemperaturesSecondToLastWeek[["2012-05-21"]] temperaturesMay21reproj <- raster(x = extent(testDataValidation@sp), nrows = nrow(temperaturesMay21), ncols = ncol(temperaturesMay21), crs = crs(indiaPolygonsOtherCRS)) temperaturesMay21reproj <- raster::projectRaster(from = temperaturesMay21, to = temperaturesMay21reproj) recordedTemperaturesInMissingZone <- raster::extract(x = temperaturesMay21reproj, y = testDataValidation@sp@coords) differences <- indiaAnalysisValidation$predMoments$Mean - recordedTemperaturesInMissingZone MSPE <- mean(differences^2, na.rm = TRUE) MedSPE <- median(differences^2, na.rm = TRUE) quantile(abs(differences), probs = 0.9) # What is the coordinate of the tile that produced the largest absolute difference? testDataSpReproj <- spTransform(testDataValidation@sp, CRSobj = crs(indiaPolygons)) testDataSpReproj@coords[which.max(abs(differences)),] spObjectReprojected <- sp::spTransform(testDataValidation@sp, CRSobj = crs(indiaPolygons)) spObjectReprojected@coords[which.min(indiaAnalysisValidation$predMoments$Mean), ] fieldValues <- differences[!is.na(differences)] sqDiffRaster <- raster::rasterize(x = spObjectReprojected, y = emptyRaster, field = fieldValues^2) diffRaster <- raster::rasterize(x = spObjectReprojected, y = emptyRaster, field = fieldValues)
The previous results were all for IS-MRA. The next chunk evaluates the predictive performance of INLA-SPDE. We re-create the mesh and transformation matrix used when we fitted the model.
library(INLA) trainingDataValidation@sp <- sp::spTransform(trainingDataValidation@sp, CRSobj = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) testDataValidation@sp <- sp::spTransform(testDataValidation@sp, CRSobj = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) timeVecTraining <- (as.numeric(time(trainingDataValidation)) - min(as.numeric(time(trainingDataValidation))))/(3600*24) + 1 timeVecTest <- (as.numeric(time(testDataValidation)) - min(as.numeric(time(trainingDataValidation))))/(3600*24) + 1 knots <- seq(1, max(timeVecTraining), length = max(timeVecTraining)) mesh1 <- inla.mesh.1d(loc = knots, degree = 2, boundary = "free") ## generate space mesh mesh2 <- inla.mesh.2d(loc = trainingDataValidation@sp@coords[timeVecTraining == 1, ], cutoff = 0.01, offset = c(0.1, 0.2), max.n = 2000) # range0 and sigma0 control the prior means for the range and scale parameters. # See Lindgren INLA tutorial page 5. d <- 1 alpha <- 2 kappa <- 1 # = 1/(range parameter in my model) spatialSmoothness <- alpha - d/2 # cf p.3 INLA tutorial loghyperparaSDinMyModel <- log(10) # range0 and sigma0 seem to be the prior means... range0 <- sqrt(8 * spatialSmoothness)/kappa # sqrt(8 * spatial smoothness) / Kappa. In my model, I use 1 as prior mean for spatial range and fix smoothness at 1.5. This means Kappa = 1. sigma0 <- 1 lkappa0 <- log(8 * spatialSmoothness)/2 - log(range0) ltau0 <- 0.5*log(gamma(spatialSmoothness)/(gamma(alpha)*(4*pi)^(d/2))) - log(sigma0) - spatialSmoothness * lkappa0 ## build the spatial spde spde <- inla.spde2.matern(mesh2, B.tau = matrix(c(ltau0, -1, spatialSmoothness), 1, 3), B.kappa = matrix(c(lkappa0, 0, -1), 1,3), theta.prior.mean = c(0,0), theta.prior.prec = c(1/loghyperparaSDinMyModel^2, 1/loghyperparaSDinMyModel^2)) ## build the space time indices STindex <- inla.spde.make.index("space", n.spde = spde$n.spde, n.group = mesh1$m) ## Link data and process Atraining <- inla.spde.make.A(mesh2, loc = trainingDataValidation@sp@coords, group = timeVecTraining, group.mesh = mesh1) Atest <- inla.spde.make.A(mesh2, loc = testDataValidation@sp@coords, group = timeVecTest, group.mesh = mesh1) stackTraining <- inla.stack(data = list(y = trainingDataValidation@data$y), A = list(Atraining, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), effects = list( c(STindex, list(intercept = 1)), list(landCover3 = trainingDataValidation@data$landCover2), # The mismatch between the element name and the name of the extracted column is to reflect an error in the original fitting, where landCover2 was encoded as landCover3. list(landCover4 = trainingDataValidation@data$landCover4), list(landCover5 = trainingDataValidation@data$landCover5), list(landCover8 = trainingDataValidation@data$landCover8), list(landCover9 = trainingDataValidation@data$landCover9), list(landCover10 = trainingDataValidation@data$landCover10), list(landCover11 = trainingDataValidation@data$landCover11), list(landCover12 = trainingDataValidation@data$landCover12), list(landCover13 = trainingDataValidation@data$landCover13), list(landCover14 = trainingDataValidation@data$landCover14), list(landCover15 = trainingDataValidation@data$landCover15), list(elevation = trainingDataValidation@data$elevation), list(Aqua = trainingDataValidation@data$Aqua), list(time2 = trainingDataValidation@data$time2), list(time3 = trainingDataValidation@data$time3), list(time4 = trainingDataValidation@data$time4), list(time5 = trainingDataValidation@data$time5), list(time6 = trainingDataValidation@data$time6), list(time7 = trainingDataValidation@data$time7) ), tag = "est") stackTest <- inla.stack( data = list(y = NA), A = list(Atest, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 , 1, 1, 1, 1), effects = list( c(STindex, list(intercept = 1)), list(landCover3 = testDataValidation@data$landCover2), # The mismatch between the element name and the name of the extracted column is to reflect an error in the original fitting, where landCover2 was encoded as landCover3. list(landCover4 = testDataValidation@data$landCover4), list(landCover5 = testDataValidation@data$landCover5), list(landCover8 = testDataValidation@data$landCover8), list(landCover9 = testDataValidation@data$landCover9), list(landCover10 = testDataValidation@data$landCover10), list(landCover11 = testDataValidation@data$landCover11), list(landCover12 = testDataValidation@data$landCover12), list(landCover13 = testDataValidation@data$landCover13), list(landCover14 = testDataValidation@data$landCover14), list(landCover15 = testDataValidation@data$landCover15), list(elevation = testDataValidation@data$elevation), list(Aqua = testDataValidation@data$Aqua), list(time2 = testDataValidation@data$time2), list(time3 = testDataValidation@data$time3), list(time4 = testDataValidation@data$time4), list(time5 = testDataValidation@data$time5), list(time6 = testDataValidation@data$time6), list(time7 = testDataValidation@data$time7) ), tag = 'predictions') combinedStack <- inla.stack(stackTraining, stackTest) stackIndex <- inla.stack.index(combinedStack, "predictions")$data preds <- SPDEresult$summary.linear.predictor differencesSPDE <- preds$mean[stackIndex] - recordedTemperaturesInMissingZone MSPE_SPDE <- mean(differencesSPDE^2, na.rm = TRUE) MedSPE_SPDE <- median(differencesSPDE^2, na.rm = TRUE) smallRasterExtent <- raster::extent(c(xmin = min(smallWesternMahaPolygonEdges[ , 1]), xmax = max(smallWesternMahaPolygonEdges[ , 1]), ymin = min(smallWesternMahaPolygonEdges[ , 2]), ymax = max(smallWesternMahaPolygonEdges[ , 2]))) fieldValuesSPDE <- differencesSPDE[!is.na(differencesSPDE)] spObjectReprojectedSPDE <- sp::spTransform(testDataValidation@sp[!is.na(differencesSPDE)], CRSobj = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) emptyRaster <- raster(x = smallRasterExtent, nrows = 120, ncols = 120, crs = sp::CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) diffRasterSPDE <- raster::rasterize(x = spObjectReprojectedSPDE, y = emptyRaster, field = fieldValuesSPDE) ecolSPDE <- mapmisc::colourScale(values(diffRasterSPDE), col = "Spectral", dec = 0, breaks = 10, rev = TRUE, style = "equal", opacity = 0.5) ecolSPDE <- ecolSPDE[c("col", "breaks")] ecolINLAMRA <- mapmisc::colourScale(values(diffRaster), col = "Spectral", dec = 0, breaks = 10, rev = TRUE, style = "equal", opacity = 0.5) ecolINLAMRA <- ecolINLAMRA[c("col", "breaks")] newBoundaries <- range(c(ecolINLAMRA$breaks, ecolSPDE$breaks)) combinedBreaks <- seq(from = newBoundaries[[1]], to = newBoundaries[[2]], length.out = length(ecolINLAMRA$breaks)) jpeg(file = "outputFiles/predErrorOnMay21withTimeCovarSPDE.jpg", width = 1200, height = 1200) raster::plot(diffRasterSPDE, legend = FALSE, cex.axis = 2.5, col = ecolINLAMRA$col, breaks = combinedBreaks) plot(indiaPolygons, add = TRUE) # plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE) # text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = citiesNames, offset = 2, cex = 5, pos = 4) mapmisc::legendBreaks("bottomleft", list(breaks = round(combinedBreaks, digits = 1), col = ecolINLAMRA$col), title = NULL, bg = "white", cex = 2.5, bty = "n") mapmisc::scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3) dev.off() jpeg(file = "outputFiles/predErrorOnMay21withTimeCovar.jpg", width = 1200, height = 1200) raster::plot(diffRaster, legend = FALSE, cex.axis = 2.5, col = ecolINLAMRA$col, breaks = combinedBreaks) plot(indiaPolygons, add = TRUE) mapmisc::legendBreaks("bottomleft", list(breaks = round(combinedBreaks, digits = 1), col = ecolINLAMRA$col), title = NULL, bg = "white", cex = 2.5, bty = "n") # raster::plot(diffRaster, legend.only = TRUE, legend.width = 4, axis.args = list(cex.axis = 3)) # plot(cityPoints, pch = 15, col = "red", cex = 5, add = TRUE) # text(x = cityPoints@coords[,1], y = cityPoints@coords[,2], labels = replace(citiesNames, 1, "Mumbai"), offset = 2, cex = 5, pos = 4) mapmisc::scaleBar(crs = crs(indiaPolygons), pos = "topleft", cex = 3, pt.cex = 2.2, title.cex = 3.5, seg.len = 3) dev.off()
The last graph we produced presents prediction errors for IS-MRA across the entire prediction region. We could not produce it before because we wanted to use the same colour scale as for the analogous graph for SPDE, to make them easily comparable when read side-by-side. We now move on to the main analysis, on the large dataset.
spAndTimeCorr <- lapply(list(space = list(paraName = "space.rho", smoothness = 1.5, distances = c(close = 1, far = 10)), time = list(paraName = "time.rho", smoothness = 0.5, distances = c(close = 1, far = 7))), function(component) { sapply(component$distances, function(distanceMeasure) { MRAinla::maternCov(d = distanceMeasure, rho = exp(indiaAnalysisMain$hyperMarginalMoments[component$paraName, "Mean"]), smoothness = component$smoothness, scale = 1) }) }) ######### Preparing parameter and hyperparameter moments table ######## hyperMoments <- subset(indiaAnalysisMain$hyperMarginalMoments[c("space.rho", "time.rho", "scale"), ], select = -Skewness) ## Delta method to restore the original scale hyperMoments$StdDev <- hyperMoments$StdDev * exp(hyperMoments$Mean) # mean(rho) = exp(mean(log rho) + 0.5 var(log rho)) hyperMoments$Mean <- exp(hyperMoments$Mean + 0.5 * hyperMoments$StdDev^2) hyperMoments[ , "CredInt_2.5%"] <- exp(hyperMoments[ , "CredInt_2.5%"]) hyperMoments[ , "CredInt_97.5%"] <- exp(hyperMoments[ , "CredInt_97.5%"]) covariateMoments <- indiaAnalysisMain$FEmarginalMoments combinedMoments <- rbind(covariateMoments, hyperMoments) landCoverNames <- c("Water", "Evergreen Needleleaf Forests", "Evergreen Broadleaf Forests", "Deciduous Needleleaf Forests", "Deciduous Broadleaf Forests", "Mixed Forests", "Closed Shrublands", "Open Shrublands", "Woody Savannas", "Savannas", "Grasslands", "Permanent Wetlands", "Croplands", "Urban and Built-up Lands", "Cropland/Natural Vegetation Mosaics", "Non-Vegetated Lands", "Unclassified") names(landCoverNames) <- paste("landCover", c(0:15, 255), sep = "") landCoverPositions <- grep(pattern = "landCover", x = rownames(combinedMoments)) landCoverValues <- grep(pattern = "landCover", x = rownames(combinedMoments), value = TRUE) rownames(combinedMoments) <- replace( rownames(combinedMoments), landCoverPositions, unname(landCoverNames[landCoverValues]) ) rownames(combinedMoments) <- replace( rownames(combinedMoments), match("elevation", rownames(combinedMoments)), "Elevation" ) rownames(combinedMoments) <- replace( rownames(combinedMoments), match("Aqua", rownames(combinedMoments)), "Satellite: Aqua" ) datesCovered <- paste("May", 26:31) names(datesCovered) <- paste("time", 2:7, sep = "") rownames(combinedMoments) <- replace( rownames(combinedMoments), grep(pattern = "time[0-9]$", x = rownames(combinedMoments)), unname(datesCovered[grep(pattern = "time[0-9]$", x = rownames(combinedMoments), value = TRUE)]) ) reorderedCombined <- combinedMoments[c("Intercept", "Elevation", "Satellite: Aqua", intersect(rownames(combinedMoments), landCoverNames), intersect(rownames(combinedMoments), datesCovered), rownames(hyperMoments)), ] latexCode <- Hmisc::latex(Hmisc::format.df(reorderedCombined, dec = 3), rgroup = c("", "Land cover", "Time", "Hyperparameters"), n.rgroup = c(3, 14, 6, 3), file = "outputFiles/JASAtable1.tex") ######## END: Preparing parameter and hyperparameter moments table ######## table((indiaAnalysisMain$predMoments$Mean - min(trainingDataMain@data[, "y"])) < -1) extremeValueIndices <- which(indiaAnalysisMain$predMoments$Mean < (min(trainingDataMain@data[, "y"]) - 1)) mostExtremeIndex <- which.min(indiaAnalysisMain$predMoments$Mean) predictionDataMain@data[mostExtremeIndex, ] predictionOrder <- order(indiaAnalysisMain$predMoments$Mean) predictionDataMain@sp@coords[predictionOrder[1:2],] extremeDataset <- predictionDataMain[extremeValueIndices] extremeDataset@data <- data.frame(Temperature = indiaAnalysisMain$predMoments$Mean[extremeValueIndices]) # jpeg("outputFiles/extremeValuePositions.jpeg", width = 1000, height = 1000) # stplot(extremeDataset) # dev.off() ######## Plotting predictions ######## indiaPolygons <- raster::getData(country = "IND", level = 2) trainingDataMainReproject <- trainingDataMain trainingDataMainReproject@sp <- sp::spTransform(trainingDataMain@sp, crs(indiaPolygons)) predictionDataMainReproject <- predictionDataMain predictionDataMainReproject@sp <- sp::spTransform(predictionDataMain@sp, crs(indiaPolygons)) # Predicted values outside the range of observed values will be omitted to get a more appealing visualisation. The following lines compute the number of omitted values: table(indiaAnalysisMain$predMoments$Mean[which(indiaAnalysisMain$predData$time == as.POSIXct("2012-05-28"))] <= min(indiaAnalysisMain$data$spObject@data$y[which(indiaAnalysisMain$data$time == as.POSIXct("2012-05-28"))])) table(indiaAnalysisMain$predMoments$Mean[which(indiaAnalysisMain$predData$time == as.POSIXct("2012-05-28"))] >= max(indiaAnalysisMain$data$spObject@data$y[which(indiaAnalysisMain$data$time == as.POSIXct("2012-05-28"))])) plot.INLAMRA(indiaAnalysisMain, filename = "outputFiles/predictionsIndiaDataMay28mainAnalysisJoint_Larger.jpg", type = "joint", polygonsToOverlay = indiaPolygons, control = list( graphicsEngine = jpeg, controlForScaleBar = list( pos = "topleft", cex = 2, pt.cex = 1.5, seg.len = 3 ), controlForRasterLegend = list( pos = "bottomleft", title = NULL, bg = "white", cex = 2.5 ), controlForRasterColourScale = list( col = "Spectral", breaks = seq(from = 20, to = 60, by = 5), style='fixed', rev = TRUE, dec = 1, opacity = 0.5 ), controlForRasterPlot = list( cex.axis = 2, cex.main = 3, main = NULL ), resolutionInMeters = 1000, trim = 2, timesToPlot = unique(time(predictionDataMain)) ), width = 1100, height = 1200 ) # There was a bug in the estimation of the SDs in an early version of the software: to get the real value, we must apply the following transformation: # sqrt(indiaAnalysisMain$predMoments$SD - errorSD^2) + errorSD errorSD <- 0.5 # cf. vignette "ISMRAmainAnalysis_Script.Rmd" indiaAnalysisMainCopy <- indiaAnalysisMain indiaAnalysisMainCopy$predMoments$SD <- sqrt(indiaAnalysisMain$predMoments$SD - errorSD^2) + errorSD plot.INLAMRA(indiaAnalysisMainCopy, filename = "outputFiles/predictionsIndiaDataMay28mainAnalysisSDs_LargerFixed.jpg", type = "SD", polygonsToOverlay = indiaPolygons, control = list( graphicsEngine = jpeg, controlForScaleBar = list( pos = "topleft", cex = 2, pt.cex = 1.5, seg.len = 3 ), controlForRasterLegend = list( pos = "bottomleft", title = NULL, bg = "white", cex = 2.5 ), controlForRasterColourScale = list( col = "Spectral", breaks = c(0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5), style='fixed', rev = TRUE, dec = 1, opacity = 1 ), controlForRasterPlot = list( cex.axis = 2, cex.main = 3, main = NULL ), resolutionInMeters = 1000, timesToPlot = unique(time(predictionDataMain)) ), width = 1100, height = 1200 ) plot.INLAMRA(indiaAnalysisMain, filename = "outputFiles/predictionsIndiaDataMay28mainAnalysisTraining_Larger.jpg", type = "training", polygonsToOverlay = indiaPolygons, control = list( graphicsEngine = jpeg, controlForScaleBar = list( pos = "topleft", cex = 2, pt.cex = 1.5, seg.len = 3 ), controlForRasterLegend = list( pos = "bottomleft", title = NULL, bg = "white", cex = 2.5 ), controlForRasterColourScale = list( col = "Spectral", breaks = seq(from = 20, to = 60, by = 5), style='fixed', rev = TRUE, dec = 0, opacity = 0.5 ), controlForRasterPlot = list( cex.axis = 2, cex.main = 3, main = NULL ), resolutionInMeters = 1000, timesToPlot = unique(time(predictionDataMain)) ), width = 1100, height = 1200 )
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.