--- title: 'Exercise: Regional mapping of climate variables (1) Data exploration' author: "D G Rossiter" date: "`r Sys.Date()`" output: html_document: toc: TRUE toc_float: TRUE theme: "lumen" number_section: FALSE fig_height: 4 fig_width: 6 fig_align: 'center' --- ```{r global_options1, include=FALSE} knitr::opts_knit$set(root.dir = '/Users/rossiter/data/edu/dgeostats/ex/ds/中国天气') writeLines(capture.output(sessionInfo()), "sessionInfo.txt") ``` ```{r global_options2, include=FALSE} knitr::opts_chunk$set(fig.width=4, fig.height=6, fig.align='center', fig.path='kgraph/', warning=FALSE, message=FALSE) # rprojroot::find_rstudio_root_file() ``` This exercise should be loaded into RStudio, in the same directory with the dataset `zhne_stations.RData`. Or, you can load this from another directory and edit the paths to the dataset. # Task 0 -- packages Load the `rgdal`, `sp`, `gstat` and `ggplot2` packages. During the following tasks, load additonal packages as needed. ```{r} library(sp); library(rgdal); library(gstat); library(ggplot2) ``` # Task 1 -- Dataset The long-term climate dataset contains 30-year averages [1981-2010] of (1) temperature annual average, minimum, maximum, temperature extreme minimum and maximum; (2) precipitation annual average, minimum, maximum at a set of stations in north-central and northeast China. Data were downloaded by 盛美玲 from the National Meteorological Information Center ^[http://data.cma.cn/site/index.html] as comma-separated (CSV) files and prepared for R by the instructor. A vector GIS coverage of administrative boundaries is also included in the R dataset `zhne_climate.RData`. The instructor prepared a 4x4 km horizontal resolution grid covering the area to be mapped. This has the elevation (m.a.s.l 海拔) according to the World 1 km SRTM improved by CGIAR ^[http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1]. Additional covariates that might help explain climate were added, as explained below: population density and a terrain index. ## Loading the dataset Load the climate stations dataset and the administrative boundaries, listing the R names of the loaded objects: ```{r load.points} load("./zhne_climate.RData", verbose=TRUE) ``` The loaded files are: * `zhne` -- observation points in WGS84 long/lat CRS, with attributes and covariates * `zhne.m` -- observation points in a metric projection (Lambert Equal Area) also on WGS84 * `zhne.df` -- same, as a data frame, with coördinates as data fields * `zh.crs` -- the Coördinate Reference System (CRS) of the `*.m` objects * `zhne.adm1.m` -- vector file of first-level administrative units covering the study area, in the Lambert CRS * `zhne.grid.m` -- 4 x 4 km resolution grid covering the area to be mapped, in the Lambert CRS * `zhne.grid.df` -- same, as a data frame, with coördinates as data fields ## Fixing the encoding This `.RData` file was saved in the UTF-8 encoding. If your system uses another encoding (typical for Chinese Windows operating systems), you have to inform R of the original encoding. It will then re-encode character strings to your system's encoding. Here is a function to loop through the fields of a dataframe and fixes them, using the `Encoding` function. I've set it up with the default orginal encoding `UTF-8` but this can be changed. For fields that are factors, it first converts them to characters, re-encodes the names, and then converts back to factors. Arguments: 1. `df` : the data frame to be re-encoded 2. `inputEncoding`, default `"utf-8"` ```{r fix.encoding} fix.encoding <- function(df, inputEncoding = "UTF-8") { numCols <- ncol(df) for (col in 1:numCols) if(class(df[, col]) == "character"){ Encoding(df[, col]) <- inputEncoding } else if(class(df[,col]) == "factor") { tmp <- as.character(df[,col]) Encoding(tmp) <- inputEncoding df[,col] <- as.factor(tmp) } return(df) } ``` Apply this to the data frames. ```{r} zhne.df <- fix.encoding(zhne.df) # only the dataframe part of the Spatial* objects zhne.m@data <- fix.encoding(zhne.m@data) zhne.adm1.m@data <- fix.encoding(zhne.adm1.m@data) ``` Check that the re-encoding worked. The Chinese characters should display correctly. ```{r} head(zhne.df$station.name) levels(zhne.df$province) levels(zhne.adm1.m$NAME_HZ) ``` ## Climate variables The climate variables are all 30-year summaries: * temperature annual average: `t.avg` * temperature annual minimum: `t.min` * temperature annual maximum: `t.max` * temperature extreme minimum: `t.ext.min` * temperature extreme maximum: `t.ext.max` * precipitation annual average: `p.avg` * precipitation annual minimum: `p.avg.min` * precipitation annual maximum: `p.avg.max` ## Covariables In addition there are covariables that can be used in mapping: * pop15: population density, average in 15' surrounding grid * pop2pt5: population density, average in 2.5' surrounding grid * mrvbf: multiresolution index of valley bottom flatness (MRVBF) See details on these next, they will not be used until the mapping part of the exercise. ### 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 1981-2010. 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 ### 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, J. C., & Dowling, T. I. (2003). A multiresolution index of valley bottom flatness for mapping depositional areas. Water Resources Research, 39(12), ESG4-1-- ESG4-13. https://doi.org/10.1029/2002WR001426] identifies valley bottoms based on their topographic signature as flat low-lying areas, at increasingly-broad scales, and combines these into a single index. This was computed in QGIS by the instructor from the 4km DEM, using the SAGA GIS implementation of the MRVBF algorithm. One parameter for the MRVBF calculation 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). ## Geographic extent Display the bounding box and CRS of the climate stations, in both CRS. ```{r} bbox(zhne) bbox(zhne.m) proj4string(zhne) proj4string(zhne.m) ``` Q: Why are the numbers in the bounding box different for `zhne` and `zhne.m`? Q: What is the centre of the Lambert Equal Area projection? Why was this chosen? Display the bounding box and CRS of the grid: ```{r} bbox(zhne.grid.m) proj4string(zhne.grid.m) ``` Q: Which of the point datasets does this match? # Task 2 -- Dataset summary ## Data structure Display the structure of the observation points in a metric projection. ```{r} str(zhne.m) ``` Q: Which slot (marked with `@`) contains the attribute data? Q: Which slot contains the coördinates? Q: How many observations are there? Q: How many attributes at each observation point? ## Summary Display a summary of the climate attributes: minimum, first quartile, median, mean, third quartile, maximum. ```{r} summary(zhne.df) ``` Q: How many stations are in Shandong? Q: What is the elevation range of the stations? Q: What is the median 30-year average annual temperature of the stations? Find the stations with 30-year average annual temperatures below 0C and display their names and locations: ```{r} table(zhne.df$t.avg < 0) zhne.df[(zhne.df$t.avg < 0), c("province", "station.name", "long", "lat", "elev.m", "t.avg")] ``` Q: Which administrative area has the most cold stations? Find the coldest station: ```{r} (ix <- which.min(zhne.df$t.avg)) zhne.df[ix, c("province", "station.name", "long", "lat", "elev.m", "t.avg")] ``` Q: Is this the northern-most station? The highest elevation station? # Task 3 -- Geographic locations Display the locations of the climate stations, within the administrative boundaries, using the metric coördinates. ```{r locations,fig.width=12, fig.height=12} plot(coordinates(zhne.m), asp=1, xlab="E", ylab="N", pch=20, cex=0.6, col=as.numeric(zhne.m$province.en) ) text(coordinates(zhne.m), labels=zhne.m$station.id, cex=0.6, pos=4) plot(zhne.adm1.m, add=TRUE, border="darkgray", lwd=1.2) grid() ``` Q: Which first-level administrative division has the densest network of observations? Find the northern-most station and display its attributes: ```{r} ix <- which.max(zhne.df$N) zhne.df[ix,] ``` Q: What is the name of this station? Do you know anything about this location? Q: Which first-level administrative divisions have no observations in this dataset, but are shown in the adminstrative units map? (We will use these to see how each model extrapolates.) # Task 4 -- Postplot of a climate variable In the in-class exercise we use `t.avg`. Display a postplot of its values at each station, size proportional to the value, coloured by elevation. ```{r postplot,fig.width=10, fig.height=10} summary(zhne.df$t.avg) require(ggplot2) ggplot(data=zhne.df) + aes(x=E, y=N) + geom_point(aes(size=t.avg, colour=elev.m), shape=20) + xlab("E") + ylab("N") + coord_fixed() ``` # Task 5 -- Relation of target and predictors ## Relation with N, E, elevation Display scatterplots of the selected climate attribute against possible predictors N, E, elevation. ```{r scatterplots,fig.width=8, fig.height=6} p1 <- ggplot() + geom_point(aes(x=elev.m, y=t.avg, colour=province.en), data=zhne.df ) p2 <- ggplot() + geom_point(aes(x=E, y=t.avg, colour=province.en), data=zhne.df ) + xlab("Easting") p3 <- ggplot() + geom_point(aes(x=N, y=t.avg, colour=province.en), data=zhne.df) + xlab("Northing") print(p3) print(p2) print(p1) ``` Q: Describe the relations between the climate attributes and the predictor. -- Do they appear to have a simple functional form, e.g., linear? -- Do they appear to be the same in all provinces? Q: Do any points not fit the overall pattern for one or more of the predictors? If so, write R code to identify it. Does this seem to be a correct point, i.e., not a database error? Why or why not? ## Relation with covariates Display scatterplots of the selected climate attribute against possible predictors population density and MRVBF. ```{r scatterplots.covariates,fig.width=8, fig.height=6} p1 <- ggplot() + geom_point(aes(x=pop15, y=t.avg, colour=province.en), data=zhne.df ) + xlab("Population density (15' grid)") p2 <- ggplot() + geom_point(aes(x=pop2pt5, y=t.avg, colour=province.en), data=zhne.df ) + xlab("Population density (2.5' grid)") p3 <- ggplot() + geom_point(aes(x=mrvbf, y=t.avg, colour=province.en), data=zhne.df) + xlab("MRVBF") print(p1) print(p2) print(p3) ``` Q: Describe the relations between the climate attributes and the predictor. -- Do they appear to have a simple functional form, e.g., linear? -- Do they appear to be the same in all provinces? # Task 6 -- Inter-relation of predictors Display pairwise scatterplots of the possible predictors. ```{r parwise,fig.width=12, fig.height=12} names(zhne.df) preds <- c(17,16,9, 18:20 ) pairs(zhne.df[, preds], col=zhne.df$province.en, pch=20, cex=0.6) ``` Q: Which pairs of predictors are closely related? Somewhat? Not at all? # Task 7 -- Grid layers ## Elevation Display the elevation as a map. ```{r show.DEM, fig.width=10, fig.height=10} ggplot(data=zhne.grid.m.df) + geom_point(aes(E, N, color=elev.m)) + xlab("E") + ylab("N") + coord_fixed() + scale_colour_distiller(name="m.a.s.l", space="Lab", palette="YlGnBu") ``` Q: Describe the main features of the relief. ## Population density Display the population density. ```{r show.popd, fig.width=12, fig.height=6} 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) ``` Q: Which map shows the highest densities? Why? Q: Which map shows more detail? Why? ## MRVBF Display the MRVBF: ```{r show.mrvbf, fig.width=8, fig.height=8} p1 <- spplot(zhne.grid.m, zcol="mrvbf", main="MRVBF", pretty=TRUE, at=seq(0, 6, length=16), col.regions=topo.colors(16)) print(p1) ``` Q: Which areas are shown as the flattest? Which with the most relief?