--- title: 'Assignment: Regional mapping of climate variables: setting up the dataset' author: "D G Rossiter" date: '`r Sys.Date()`' output: html_document: code_folding: show fig_align: center fig_height: 4 fig_width: 6 number_section: yes theme: lumen toc: yes toc_float: yes pdf_document: toc: yes word_document: toc: yes --- ```{r global_options, include=FALSE} knitr::opts_chunk$set(fig.width=4, fig.height=6, fig.align='center', fig.path='kgraph/', warning=FALSE, message=FALSE) writeLines(capture.output(sessionInfo()), "sessionInfo.txt") ``` This script prepares long-term climate records from a set of stations in China for analysis with the methods of the tutorial "Applied geostatistics: Regional mapping of climate variables from point samples", which uses a dataset from the northeastern USA. # Dataset source This dataset contains 30-year averages of (1) temperature annual average, minimum, maximum, temperature extreme minimum and maximum; (2) precipitation annual average, minimum, maximum. Data were downloaded from the National Meteorological Information Center ^[http://data.cma.cn/site/index.html] as comma-separated (CSV) files by 盛美玲. These files are: * `temperature_stn_1981_2010.csv`: temperature * `precip_stn_1981_2010.csv` : precipitation # Preparation Packages: ```{r packages} library(sp) # spatial data structure for R library(rgdal) # coordinate reference systems, interface to geographic data library(rgeos) # manipulate geometry library(dplyr) # modern data wrangling library(stringr) # modern string manipulation ``` # Import of the temperature dataset. The station and province names use Chinese characters, so we need to identify the encoding. Here are the encoding codes: ```{r encodings} getOption("encoding") # GB-series codes 国家标准 iconvlist()[ix <- grep("GB", iconvlist(), fixed=TRUE)] # UTF-series code: Unicode Transformation Format iconvlist()[ix <- grep("UTF", iconvlist(), fixed=TRUE)] ``` We guess that the encoding is the common `GB18030`; if this does not work, we would try others, until the characters are legible. Read the temperature dataset; guess the encoding: ```{r read.temp} zhne <- read.csv("./temperature_stn_1981_2010.csv", fileEncoding="GB18030", stringsAsFactors=FALSE) str(zhne) ``` Note: we use the R name `zhne` = 中国东北 because most of these stations are in the northeast of China (although the area is wider than the common definition of 中国东北). Recode the province names to Latin alphabet, keep the originals: ```{r provinces} zhne$province <- as.factor(zhne$province) levels(zhne$province) # northeast PRC zhne$province.en <- recode(zhne$province, "北京"="Beijing", "天津"="Tianjin", "河北"="Hebei", "山西"="Shanxi", "内蒙古"="NeiMenggu", "辽宁"="Liaoning", "吉林"="Jilin", "黑龙江"="Heilongjiang", "山东"="Shandong") table(zhne$province.en) ``` After some experiments, it became clear that the coördinates are in the format `dd.mm`, not decimal degrees. So we have to convert these to decimal degrees. ```{r} ddmm.to.dd <- function(ddmm) { .deg <- ddmm%/%1 .min <- (ddmm%%1)*100 .dd <- .min/60 return(.deg + .dd) } head(zhne$longitude) head(zhne$latitude) zhne$longitude <- ddmm.to.dd(zhne$longitude) zhne$latitude <- ddmm.to.dd(zhne$latitude) head(zhne$longitude) head(zhne$latitude) ``` Convert to a `SpatialPointsDataFrame` by specifying the coördinates, keep them also as fields in the data frame. ```{r convert.to.sp} names(zhne) coordinates(zhne) <- c("longitude", "latitude") # they were removed as fields, add them back zhne$long <- coordinates(zhne)[,1]; zhne$lat <- coordinates(zhne)[,2] str(zhne) ``` Find the time period for each station's records: ```{r station.time} table(zhne$start.year) table(zhne$end.year) ``` Use just the stations with the full 30-year record: ```{r long.term} (n <- length(ix <- which(zhne$start.year != 198101))) if (n > 0) { zhne <- zhne[-ix,] } dim(zhne) ``` Remove the date columns: ```{r remove.dates.fields.1 } c.ix <- which(!is.na(match(names(zhne), "start.year"))) zhne@data <- zhne@data[, -c.ix] c.ix <- which(!is.na(match(names(zhne), "end.year"))) zhne@data <- zhne@data[, -c.ix] c.ix <- which(!is.na(match(names(zhne), "缺测时段"))) zhne@data <- zhne@data[, -c.ix] ``` Shorten the names of the data fields, for easier typing in R models: ```{r} names(zhne) names(zhne)[1] <- "station.id" names(zhne)[3:7] <- c("t.avg", "t.avg.max", "t.avg.min", "t.ext.max", "t.ext.min") names(zhne)[10] <- "elev.m" names(zhne) ``` Assign the Coordinate Reference System (CRS) to the spatial dataset. Clearly, the coordinates are long/lat: ```{r summary.coords} bbox(zhne) ``` Since this is older data, we do not think this is WGS84 for the long/lat. 盛美玲 could not find the information about what datum/ellipse was used for the lon/lat in these datasets. We suspect that the Krasovsky 1942 ellipse is correct. Find it in the EPSG database and set the CRS accordingly, ```{r set.crs} getPROJ4VersionInfo() str(ellps <- projInfo("ellps")) (ix <- grep("Krassovsky, 1942", ellps$description)) ellps[ix,] ## this is the one we want str(proj <- projInfo("proj")) (ix <- grep("Lat/long", proj$description)) proj[ix,] # all aliases # if we think it is Krassovsky (proj4string(zhne) <- "+proj=longlat +ellps=krass") # if we think it is WGS84 # (proj4string(zhne) <- "+init=epsg:4326") bbox(zhne) ``` # Import of the precipitation dataset Read from the CSV file: ```{r read.precip} ds.p <- read.csv("./precip_stn_1981_2010.csv", fileEncoding="GB18030", stringsAsFactors=FALSE) str(ds.p) names(ds.p)[1] <- "station.id" ``` Limit to the records with a full 30-year series: ```{r} n <- length(ix <- which(ds.p$start.year != 198101)) if (n > 0) { ds.p <- ds.p[-ix,] } dim(ds.p) ``` # Add precipitation data to the temperature dataframe Do the two datasets have the same stations? ```{r} table(zhne$station.id==ds.p$station.id) table(zhne$station.name==ds.p$station.name) ``` Yes, so we can add the precipitation records to the temperature dataset, to have a single "climate" set. ```{r} names(ds.p) # note there are no 'extremes' here, as in the temperature dataset names(ds.p)[3:5] <- c("p.avg", "p.avg.max", "p.avg.min") dim(zhne@data) zhne@data <- cbind(zhne@data, ds.p[,3:5]) dim(zhne) names(zhne) ``` # Deal with duplicate locations Several stations have two records, covering two different sections of the 30-year period: ```{r zero.dist} dim(zd <- zerodist(zhne))[1] ``` There are `r dim(zd)[1]` pairs. Note that they follow each other in the list. ```{r show.zero} for (i in 1:dim(zd)[1]) { print(zhne@data[zd[i,], c(1,2,13,14)]) } ``` For each, find the number of years in each of the two records, and build a weights matrix showing the weight to give each of the two parts, when averaging the climate attributes. ```{r get.wts} wts <- matrix(0, nrow=dim(zd)[1], ncol=2) for (.i in 1:dim(zd)[1]) { # both records of the pair have the same string, use the first .t <- zhne@data[zd[.i,1], "标准值统计时段"] # two time periods .t <- strsplit(.t, "/", fixed=TRUE)[[1]] # beginning and end of each period .t1 <- strsplit(.t[1], "-", fixed=TRUE)[[1]] .t2 <- strsplit(.t[2], "-", fixed=TRUE)[[1]] # first period: beginning and end year, month .y <- as.numeric(substr(.t1, 1, 4)) .m <- as.numeric(substr(.t1, 5, 6)) # first period: number of years, months .mm <- .m[2]-.m[1]+1; .yy <- .y[2]-.y[1] # first period: total months, to use as weight .w1 <- .yy*12 + .mm # same for 2nd period .y <- as.numeric(substr(.t2, 1, 4)) .m <- as.numeric(substr(.t2, 5, 6)) .mm <- .m[2]-.m[1]+1; .yy <- .y[2]-.y[1] .w2 <- .yy*12 + .mm .sum <- .w1+.w2 wts[.i,1] <- .w1/.sum; wts[.i,2] <- .w2/.sum } ``` Now use these weights to average the two records for each pair. ```{r apply.wts} names(zhne@data) for (.i in 1:dim(zd)[1]) { .v1 <- zhne@data[zd[.i,1],] .v2 <- zhne@data[zd[.i,2],] for (.j in c(3:7, 15:17)) { .x1 <- .v1[1,.j]; .x2 <- .v2[1,.j] # overwrite values in the first pair with the weighted average zhne@data[zd[.i,1],.j] <- wts[.i,1]*.x1 + wts[.i,2]*.x2 } } ``` Remove the second record of the pair: ```{r} dim(zhne) zhne <- zhne[-zd[,2],] dim(zhne) ``` Check that no duplicates remain: ```{r} zerodist(zhne) ``` Remove the time period fields: ```{r remove.date.fields.2} c.ix <- which(!is.na(match(names(zhne), "V_TIME_AVAILA"))) zhne@data <- zhne@data[, -c.ix] c.ix <- which(!is.na(match(names(zhne), "标准值统计时段"))) zhne@data <- zhne@data[, -c.ix] names(zhne) ``` # Define a metric CRS and transform to this To build geostatistical models we need a metric CRS, i.e., with more or less true distances, which is not the case with long/lat. Define a metric CRS, for this approx. square area, we select Lambert Azimuthal Equal Area (see Snyder references for a discussion of the properties of various projections). First, look for Chinese CRS defined in the EPSG database^[https://www.epsg-registry.org]: ```{r} epsg <- make_EPSG() names(epsg) epsg[grep("China", epsg$note), 1:2] ``` These are too far south, and none are Lambert. Second, look for Lambert Azimuthal Equal Area projections: ```{r} epsg[grep("laea", epsg$prj4),1:2] ``` None for China or even Asia. So we define a custom CRS. The Lambert Azimuthal Equal Area projection is defined as follows: * +proj=laea : name of the projection in proj4 * +lon_0= : Longitude of projection center. Defaults to 0.0. * +lat_0= : Latitude of projection center. Defaults to 0.0. * +ellps= : See proj -le for a list of available ellipsoids, Defaults to “GRS80”. * +R= : Radius of the sphere given in meters. If used in conjunction with +ellps +R takes precedence. * +x_0= : False easting from projection center. Defaults to 0.0. * +y_0= :False northing from projection center. Defaults to 0.0 ## Select an ellipsoid These are possible ellipsoids: ```{r} sort(as.character(projInfo("ellps")$name)) ``` We select `WGS84`. ## Find an appropriate projection centre We use the centre of the bounding box as the projection centre: ```{r} bbox(zhne) (lon.0 <- round(mean(bbox(zhne)[1,]),1)) (lat.0 <- round(mean(bbox(zhne)[2,]),1)) ``` ## proj4 string for the coördinate reference system No reason to use false N or E. ```{r} (zh.crs <- paste0("+proj=laea +ellps=WGS84 +lon_0=", lon.0, " +lat_0=", lat.0)) ``` # Transform the point coverages Change the coördinate names to show that this is a metric, not geographic system: ```{r} zhne.m <- spTransform(zhne, zh.crs) bbox(zhne.m) coordnames(zhne.m) <- c("E", "N") ``` Show the points with average annual temperature and precipitation: ```{r show.maps, fig.width=7, fig.height=7} spplot(zhne.m, zcol="t.avg", key.space="right", pch=20, cex=0.5, main="Annual average temperature, 30-year average") spplot(zhne.m, zcol="p.avg", key.space="right", pch=20, cex=0.5, main="Annual average precipitation, 30-year average") ``` Also make a dataframe version, with the transformed coördinates as data fields, some modelling approaches only use data fields. ```{r} zhne.df <- as.data.frame(zhne.m) names(zhne.df) ``` # Export point coverages Export this dataset as a geopackage, to check the locations within a GIS. For this, use the untransformed points. Note that geopackages can contain multiple coverages. ```{r} writeOGR(zhne, dsn="./climate.gpkg", layer="climate", driver="GPKG", overwrite_layer=TRUE) ``` Difficult to check, because no reliable background map. # Political boundaries These will be used to clip the covariates, and also to display the station locations. Lots of problems with the GADM files ^[https://gadm.org/download_country_v3.html], so tried from Harvard ^[http://worldmap.harvard.edu/data/geonode:Provinces_1997]: ```{r import.prov} zh.tm <- readOGR(dsn="./Provinces_1997", layer="Provinces_1997") summary(zh.tm) ``` This is in a custom transmercator projection. ```{r limit.prov} levels(zh.tm$NAME_HZ) # Include 陝西 and "河南", in same climate region, for extrapolation provinces.with.stations <- c("内蒙古", "北京", "吉林", "天津", "山東", "山西", "河北", "遼寧", "黑龍江", "陝西", "河南") zhne.adm1.m <- zh.tm[zh.tm$NAME_HZ %in% provinces.with.stations, ] bbox(zhne.adm1.m) # drop unused levels zhne.adm1.m@data$NAME_HZ <- factor(zhne.adm1.m@data$NAME_HZ) zhne.adm1.m@data$NAME_PY <- factor(zhne.adm1.m@data$NAME_PY) levels(zhne.adm1.m$NAME_HZ) levels(zhne.adm1.m$NAME_PY) ``` See how these are in geographic coords: ```{r transform.ne.polys.sp} tmp <- spTransform(zhne.adm1.m, CRS("+proj=longlat +datum=WGS84")) bbox(tmp) ``` This agrees with a visual estimate of the northmost point of China, on Google Earth. Transform to the CRS used in this project and give a more useful name: ```{r transform.ne.polys, fig.width=8, fig.height=6} zhne.adm1.m <- spTransform(zhne.adm1.m, zh.crs) bbox(zhne.adm1.m) spplot(zhne.adm1.m, zcol="NAME_PY") ``` Some stations are not in their declared boundaries: ```{r show.noprov.stations, fig.height=9, fig.width=9} ix <- over(zhne.m, zhne.adm1.m) length(no.province <- which(is.na(ix$NAME_PY))) zhne.m@data[no.province, c("station.id", "province", "station.name", "long", "lat")] plot(coordinates(zhne.m), asp=1, xlab="E", ylab="N", pch=20, cex=0.6, col=ifelse(is.na(ix$NAME_PY), "red", "green")) text(coordinates(zhne.m), labels=ifelse(is.na(ix$NAME_PY), zhne.m$station.id, ""), cex=0.6, pos=4) plot(zhne.adm1.m, add=TRUE, border="darkgray", lwd=1.2) grid() ``` These are just outside the simplified boundaries. But on overlay with a coarse 4 km grid they should match a cell and get attribute values, so keep them. # Save stations and provinces as R objects Save as an R object: ```{r save.r} save(zhne, zhne.m, zhne.df, zhne.adm1.m, zh.crs, file="./zhne_stations.RData") ``` # References * Bugayevskiy, L. M., & Snyder, J. P. (1995). Map projections: a reference manual. Taylor & Francis. * Snyder, J. P. (1987). Map projections: a working manual. Retrieved from https://pubs.er.usgs.gov/publication/pp1395 * Snyder, J. P., & Voxland, P. M. (1989). An album of map projections. Retrieved from https://pubs.er.usgs.gov/publication/pp1453