--- title: "Tutorial: Regional mapping of climate variables from point samples (Northeast China) -- setting up additional covariates" author: "D G Rossiter" date: "`r Sys.Date()`" output: html_document: toc: TRUE toc_float: TRUE theme: "lumen" code_folding: show number_section: TRUE fig_height: 4 fig_width: 6 fig_align: 'center' bibliography: /Users/rossiter/data/edu/dgeostats/ex/ex.bib --- ```{r global_options, include=FALSE} knitr::opts_chunk$set(fig.width=4, fig.height=6, fig.align='center', fig.path='kgraph/add', warning=FALSE, message=FALSE) writeLines(capture.output(sessionInfo()), "sessionInfo_zhClimateSetupAdditionalCovariates.txt") ``` # Introduction In this document we set up additional covariates (other than N, E, elevation which are already known), which can be used in mapping methods. Packages used in this tutorial: ```{r} library(raster) library(sp) library(rgdal) library(rgeos) library(ggplot2) ``` # Stations and interpolation grid These were developed in `ZhClimateAssignment_DatasetSetup.Rmd`, we use them here to (1) limit the bounding box of covariates; (2) extract values of the covariates at stations and across the grid. At the end of this tutorial we save these with the new covariates, ready for modelling. Load the weather station locations and attributes, and make a local variable with its CRS: ```{r} list.files(pattern="RData$") load(file="./zhne_stations.RData", verbose=TRUE) names(zhne.df) bbox(zhne.m) bbox(zhne) zh.crs ``` Load the DEM covering the study area and make a local variable with its extent. Also make a `SpatialPoints` version for functions that can not deal with rasters. ```{r} load(file="./dem_ne_4km.RData", verbose=TRUE) class(dem.ne.m.sp) names(dem.ne.m.sp) summary(dem.ne.m.df) ``` Copy the grids to a more obvious name, and use this when saving: ```{r rename.grid} zhne.grid.m <- dem.ne.m.sp zhne.grid.df <- dem.ne.m.df zhne.grid.ll <- spTransform(zhne.grid.m, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")) (dem.extent <- extent(zhne.grid.m)) # bounding box of the DEM, metric coords rm(dem.ne.m.sp) rm(dem.ne.m.df) ``` Convert the extent of the raster to a polygon (4-corners); assign its CRS; convert to unprojected on the WGS84 ellipse; this is then the extent with which to crop rasters in unprojected CRS to our study area. ```{r} e.p <- as(dem.extent, "SpatialPolygons") proj4string(e.p) <- zh.crs e.p.wgs84 <- spTransform(e.p, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")) # this was too small, I can not figure out why # (crop.box.84 <- bbox(e.p.wgs84)) # # so use the original lat/long from the NE China DEM (crop.box.84 <- bbox(zhne.grid.ll)) extent(crop.box.84) ``` # Population Density The concept here is that urban areas affect the local climate (typically making it warmer), and that very sparsely-populated rural areas (typically cooler, more precipitation as snow). The Socioeconomic Data and Applications Center (sedac) at Columbia University (USA) [^5] is a data centre in NASA's Earth Observing System Data and Information System (EOSDIS). One product available there is the Gridded Population of the World (GPW), v4 [^6] for several dates. We select the year 2000, since our climate data is 1971-2000. There are two grid resolutions: 2.5' (approx. 5 km) and 15' (approx. 30 km) resolutions; we try them both, because they give two scales of possible population effects. The units of measure in both cases are units: persons per $\mathrm{km}^2$. [^5]: https://sedac.ciesin.columbia.edu [^6]: https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev11 ## Download From the web page, we download the GeoTIFF files `gpw_v4_population_density_rev11_2000_15_min.tif` and `gpw_v4_population_density_rev11_2000_2pt5_min.tif` into our data directory. These are gridded dataset and can most conveniently be read and manipulated by the `raster` package. Load into the workspace: ```{r} pop2pt5.world <- raster("/Users/rossiter/data/edu/dgeostats/ex/ds/NEweather/gpw_v4_population_density_rev11_2000_2pt5_min.tif") pop15.world <- raster("/Users/rossiter/data/edu/dgeostats/ex/ds/NEweather/gpw_v4_population_density_rev11_2000_15_min.tif") proj4string(pop2pt5.world) res(pop2pt5.world) proj4string(pop15.world) res(pop15.world) ``` Both of these are in unprojected geographic coordinates. ## Restrict to the study area These maps cover the whole world; we only need our study area. Restrict the rasters to the bounding box of the study area, which we have from the DEM used in the exercise. We use an unprojected version of the box. Now crop: ```{r} pop15.ne <- raster::crop(pop15.world, crop.box.84) pop2pt5.ne <- raster::crop(pop2pt5.world, crop.box.84) names(pop15.ne) <- "pop15" names(pop2pt5.ne) <- "pop2pt5" extent(pop15.ne) extent(pop2pt5.ne) summary(pop15.ne) summary(pop2pt5.ne) ``` ## Match the projection of the DEM and stations The CRS must match the weather station and DEM CRS. ```{r} pop15.ne.m <- projectRaster(pop15.ne, crs=zh.crs, method="bilinear") nrow(pop15.ne.m); ncol(pop15.ne.m) res(pop15.ne.m) # approx. 21 x 28 km pop2pt5.ne.m <- projectRaster(pop2pt5.ne, crs=zh.crs, method="bilinear") nrow(pop2pt5.ne.m); ncol(pop2pt5.ne.m) res(pop2pt5.ne.m) # approx. 3.5 x 4.6 km ``` Check for impossible values, i.e., $< 0$, artefacts of the projection. ```{r} summary(pop15.ne.m) pop15.ne.m@data@values <- pmax(pop15.ne.m@data@values, 0, pop15.ne.m@data@values) summary(pop2pt5.ne.m) pop2pt5.ne.m@data@values <- pmax(pop2pt5.ne.m@data@values, 0, pop2pt5.ne.m@data@values) ``` ## Extract the values at the climate stations for modelling Use the `sp::over` function; this does not work on rasters, so the raster must be converted to an `sp` object. Determine the population densities of the raster cell in which each station falls, and add this as an attribute of the station. ```{r missing.pop15} bbox(zhne.m) bbox(pop15.ne.m) ix.vals <- as.vector(over( zhne.m, as(pop15.ne.m, "SpatialPixelsDataFrame"))) head(ix.vals) zhne.df$pop15 <- ix.vals$pop15 summary(zhne.df$pop15) ix.na <- which(is.na(zhne.df$pop15)) zhne.df[ix.na, c("province", "station.name", "long", "lat")] coordinates(zhne.m[ix.na,]) ``` One cell is missing a value, this is 成山头 at the extreme E of 山东. Get the population density from the adjacent cell to the west, for 15' population this is about 14 km (half of a grid cell) inland. ```{r} pt <- as.data.frame(zhne.m[ix.na,]) pt$E <- pt$E - 14000 coordinates(pt) <- ~ E+ N proj4string(pt) <- proj4string(pop15.ne.m) (val <- as.vector(over(pt, as(pop15.ne.m, "SpatialPixelsDataFrame")))) zhne.df[ix.na, "pop15"] <- val ``` ```{r missing.pop2pt5} ix.vals <- as.vector(over( zhne.m, as(pop2pt5.ne.m, "SpatialPixelsDataFrame"))) zhne.df$pop2pt5 <- ix.vals$pop2pt5 summary(zhne.df$pop2pt5) ix.na <- which(is.na(zhne.df$pop2pt5)) zhne.df[ix.na,c("province", "station.name", "long", "lat")] rm(ix.vals) ``` No missing values, good. ```{r} names(zhne.df) zhne.df[which.min(zhne.df$pop15),c("station.name", "province")] zhne.df[which.max(zhne.df$pop15),c("station.name", "province")] zhne.df[which.min(zhne.df$pop2pt5),c("station.name", "province")] zhne.df[which.max(zhne.df$pop2pt5),c("station.name", "province")] ``` There is clear evidence of dilution of density with increasing pixel size, so we have two scales of possible predictors. We see a wide spread. There are stations with almost no population nearby (e.g., remote areas of Inner Mongolia) and some with very high population density (e.g., Beijing) Copy these attributes from the dataframe to the spatial object: ```{r} zhne.m$pop15 <- zhne.df$pop15 zhne.m$pop2pt5 <- zhne.df$pop2pt5 ``` Display the relation between the target variable and log-population density: ```{r gdd.vs.pop.dens, fig.width=6, fig.height=4} plot(t.avg ~ pop15, data=zhne.df, xlab="Population density, 15'", bg=zhne.df$province, pch=21, log="x") plot(t.avg ~ pop2pt5, data=zhne.df, xlab="Population density, 2.5'", bg=zhne.df$province, pch=21, log="x") ``` There is a weak linear relation with log-population, but that is likely because of the population distribution, more people live in the warmer places. ## Transformation This predictor is highly-skewed, which will distort any method based on variances, such as regression trees. So, we transform these to logarithms, below we do the same for the prediction grid. ```{r pop.hist, fig.width=6, fig.height=4} hist(zhne.df$pop15, breaks=24) rug(zhne.df$pop15) hist(zhne.df$pop2pt5, breaks=24) rug(zhne.df$pop2pt5) ``` Add a small value to avoid taking the logarithm of 0. ```{r pop.to.log10, fig.width=6, fig.height=4} zhne.df$pop15 <- log10(zhne.df$pop15+1) zhne.df$pop2pt5 <- log10(zhne.df$pop2pt5+1) hist(zhne.df$pop15, breaks=24) rug(zhne.df$pop15) hist(zhne.df$pop2pt5, breaks=24) rug(zhne.df$pop2pt5) ``` ## Extract the values at each grid cell Add population densities to the grid. Note these are *not* the populations of the cell, rather, the density of the area surrounding and including it. ```{r} ix.vals <- over( zhne.grid.m, as(pop15.ne.m, "SpatialPixelsDataFrame")) names(ix.vals) names(zhne.grid.m) eval(parse(text=paste("zhne.grid.m$pop15 <- ix.vals$", names(ix.vals), sep=""))) names(zhne.grid.m) ix.vals <- over( zhne.grid.m, as(pop2pt5.ne.m, "SpatialPixelsDataFrame")) names(ix.vals) eval(parse(text=paste("zhne.grid.m$pop2pt5 <- ix.vals$", names(ix.vals), sep=""))) rm(ix.vals) names(zhne.grid.m) ``` Check for cells with NA's, change these to 0 (no population). ```{r} ## no NA's for 15' length(ix <- which(is.na(zhne.grid.m$pop15))) if (length(ix) > 0) zhne.grid.m[ix,"pop15"] <- 0 ## some NA's for 2.5' -- assume these are 0 pop density length(ix <- which(is.na(zhne.grid.m$pop2pt5))) if (length(ix) > 0) zhne.grid.m[ix,"pop2pt5"] <- 0 ``` Convert to logarithms: ```{r dem.pop.to.log} zhne.grid.m$pop15 <- log10(zhne.grid.m$pop15+1) zhne.grid.m$pop2pt5 <- log10(zhne.grid.m$pop2pt5+1) summary(zhne.grid.m$pop15) summary(zhne.grid.m$pop2pt5) ``` ```{r show.popdensity, fig.width=8, fig.height=4} p1 <- spplot(zhne.grid.m, zcol="pop15", main="density, log_10 per km2, 15'", pretty=TRUE, at=seq(0, 6, length=16), col.regions=topo.colors(16)) p2 <- spplot(zhne.grid.m, zcol="pop2pt5", main="density, log_10 per km2, 2.5'", pretty=TRUE, at=seq(0, 6, length=16), col.regions=topo.colors(16)) print(p1, split=c(1,1,2,1), more=T) print(p2, split=c(2,1,2,1), more=F) ``` Obvious smoothing effect at the coarser resolution. # Terrain For a good introduction to terrain analysis, see [@Wilson_Gallant_2000]. In this tutorial we have selected two terrain indices that have a relation with local topography, which is where we expect any local effects on climate. ## Preparation in QGIS Terrain covariates are all derived from the base DEM, in this example a fairly coarse $\approx 4 \times 4$ km grid. Among the programs to compute terrain covariates, SAGA GIS [^7] is one of the most comprehensive. The SAGA algorithms are also available in QGIS 3 [^8]. [^7]: http://saga-gis.org [^8]: https://www.qgis.org ## MRVBF Local terrain may influence climate. For example, narrow valleys may be "frost pockets" and often have morning ground fogs in spring and early summer. Broad plains tend to have longer-range spatial dependence than high-relief areas. The multiresolution index of valley bottom flatness (MRVBF) [@Gallant_Dowling_2003] identifies valley bottoms based on their topographic signature as flat low-lying areas, at increasingly-broad scales, and combines these into a single index. One parameter for MRVBF was adjusted from defaults: initial threshold for slope 4% (not 16%) because of the coarse resolution of the DEM. The others were default: lowness threshold 0.4, upness threshold 0.35, shape for slope 4, shape for elevation 3, maximum resolution 100%; not classified (i.e., retain the continuous value). The result was saved as a SAGA file set; these have extension `.sdat` for the data (raster), `.sgrd` for the grid parameters, `.prj` for the projection information, and `.mgrd` for the grid metadata, i.e., the processing steps. This latter is useful to confirm how the raster was created. ## Import to R We specify the full file name for the grid data, and the associated files are automatically consulted as needed: ```{r import.mrvbf, fig.height=4, fig.width=4} mrvbf <- raster("./mrvbf_ne.sdat") proj4string(mrvbf) # formatted differently but same as zh.crs extent(mrvbf) # not same as DEM image(mrvbf, asp=1, col=bpy.colors(16)) ``` These have the same CRS as the stations and DEM, but SAGA has changed the order of parameters. So, adjust the CRS to be identical: ```{r} proj4string(mrvbf) <- zh.crs ``` No need to crop, it has the same extent as the DEM, and we will extract the values at points and pixels. ## Extract the values at the climate stations for modelling Extract these and save as fields in the stations objects: Use the `sp::over` function; this does not work on rasters, so the raster must be converted to an `sp` object. ```{r} ix.vals <- over( zhne.m, as(mrvbf, "SpatialPixelsDataFrame")) names(ix.vals) names(zhne.m) eval(parse(text=paste("zhne.m$mrvbf <- ix.vals$", names(ix.vals), sep=""))) eval(parse(text=paste("zhne.df$mrvbf <- ix.vals$", names(ix.vals), sep=""))) names(zhne.m) rm(ix.vals) ``` Check for points with NA's: ```{r} length(ix <- which(is.na(zhne.m$mrvbf))) zhne.m@data[ix,] ``` Display the spatial distribution of these attributes at the stations: ```{r show.terrain.stations, fig.width=8, fig.height=4} p1 <- spplot(zhne.m, zcol="mrvbf", key.space="right", main="MRVBF") print(p1) ``` There is a good spatial distribution of these. Is there any relation with the target variable? ```{r gdd.vs.terrain, fig.width=6, fig.height=4} plot(t.avg ~ mrvbf, data=zhne.df, xlab="Multi-resolution valley bottom flatness", bg=zhne.df$province, pch=21) legend("bottomright", legend=levels(zhne.df$province.en), pch=20, col=1:4) ``` Seems very little. ## Extract the values at each grid cell Add terrain parameter to the grid. ```{r} ix.vals <- over( zhne.grid.m, as(mrvbf, "SpatialPixelsDataFrame")) names(ix.vals) eval(parse(text=paste("zhne.grid.m$mrvbf <- ix.vals$", names(ix.vals), sep=""))) rm(ix.vals) names(zhne.grid.m) ``` Check for cells with NA's: ```{r} ## no NA's for 15' length(ix <- which(is.na(zhne.grid.m$mrvbf))) ``` None, good. Display the grid: ```{r show.mrvbf, fig.width=8, fig.height=4} p1 <- spplot(zhne.grid.m, zcol="mrvbf", main="MRVBF", pretty=TRUE, at=seq(0, 6, length=16), col.regions=topo.colors(16)) print(p1) ``` # Save the extended dataset Save the stations and the interpolation grid, with the added covariates. ```{r} class(zhne.grid.m) names(zhne.grid.m) coordnames(zhne.grid.m) <- c("E", "N") zhne.grid.m.df <- as(zhne.grid.m, "data.frame") save(zhne, zhne.m, zhne.df, zh.crs, zhne.adm1.m, zhne.grid.m, zhne.grid.m.df, file="./Zh_StationsDEM_covariates.RData") ``` These are now ready for modelling and prediction. # References