---
title: "Assignment: Regional mapping of climate variables: creating a regional grid"
author: "D G Rossiter"
date: "`r Sys.Date()`"
output:
   html_document:
     theme: "lumen"
     toc: yes
     toc_float: yes
---

# World 1 km SRTM import

```{r}
library(rgdal)
library(raster)
```

The World 1 km SRTM uncompressed ASCII is 3.1 Gb. I have copied it from  http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1 to Google Drive at `https://goo.gl/6sFVwC`, then downloaded the RAR, decompressed to the `.asc`.

This seemed to be corrupted, so I got the ESRI version `SRTM_1km_GRD`.


```{r read.raster.grd}
system.time(
  dem.1km <- raster("/Users/rossiter/data/edu/dgeostats/ex/ds/SRTM_1km_GRD/srtmv4_30s"))
```

Examine the raster:

```{r}
dataType(dem.1km)
class(dem.1km)
extent(dem.1km)
projection(dem.1km)
bbox(dem.1km)   # another way to see these
proj4string(dem.1km)
res(dem.1km); nrow(dem.1km); ncol(dem.1km); ncell(dem.1km)
```


Ground resolution in km at the central latitude of our study area:

```{r}
round(res(dem.1km)[1]*(10000/90),3)
round(res(dem.1km)[1]*(10000/90)*cos(44.3*pi/180),3)
```

So about 1 km N-S, 2/3 km E-W.

Load the stations and provinces saved from `ZhClimateAssignment_DatasetSetup.Rmd`.

```{r get.stations.prov}
load("zhne_stations.Rdata", verbose=TRUE)
```

# Crop to provinces bounding box

Make a long/lat version of the study area boundary, to crop the original DEM:

```{r prov.ll}
zh.ne.ll <- spTransform(zh.ne.m, projection(dem.1km))
bbox(dem.1km); bbox(zh.ne.ll)
```

Now crop to that  bounding box:

```{r crop.prov}
system.time(
  dem.ne.ll <- crop(dem.1km, extent(zh.ne.ll))
)
names(dem.ne.ll)
names(dem.ne.ll) <- "elev.m"
# remove negative elevations, here they are artefacts of the rasterization
dem.ne.ll <- calc(dem.ne.ll, function(x) { pmax(0, x) })
summary(dem.ne.ll)
bbox(dem.ne.ll)
```

# Mask to provinces

Make a version masked to just the provinces. Keep the full bbox to compute terrain parameters at the edges of the provinces.

```{r mask.prov}
dem.ne.ll.bbox <- dem.ne.ll # save the full box
system.time(
  dem.ne.ll <- mask(dem.ne.ll, zh.ne.ll)
)
summary(dem.ne.ll)
bbox(dem.ne.ll)
image(dem.ne.ll)
```


# Project to a metric CRS

Reproject to metric projection used for climate stations, also setting the 4km resolution:

```{r reproject.dem}
system.time(
  dem.ne.m.bbox <- projectRaster(dem.ne.ll.bbox,
                          crs=proj4string(zhne.m),
                          method="bilinear",
                          res=c(4000, 4000))
)
nrow(dem.ne.m.bbox); ncol(dem.ne.m.bbox); ncell(dem.ne.m.bbox); res(dem.ne.m.bbox)
summary(dem.ne.m.bbox)
system.time(
  dem.ne.m <- projectRaster(dem.ne.ll,
                          crs=proj4string(zhne.m),
                          method="bilinear",
                          res=c(4000, 4000))
)
nrow(dem.ne.m); ncol(dem.ne.m); ncell(dem.ne.m); res(dem.ne.m)
summary(dem.ne.m)

```




# Display and save

Display the bbox, projected and aggregated DEM:

```{r display.dem.m.bbox, fig.width=8, fig.height=8}
image(dem.ne.m.bbox, col=topo.colors(24), asp=1)
```

Display the masked, projected and aggregated DEM:

```{r display.dem.m, fig.width=8, fig.height=8}
image(dem.ne.m, col=topo.colors(24), asp=1)
```

Convert to a `data.frame` for linear modelling:

```{r}
names(dem.ne.m)
names(dem.ne.m) <- "elev.m"
dem.ne.m.sp <- as(dem.ne.m, "SpatialPixelsDataFrame")
coordnames(dem.ne.m.sp) <- c("E", "N")
names(dem.ne.m.sp)
dem.ne.m.df <- as.data.frame(dem.ne.m.sp)
```

Remove the negative (impossible) elevations:

```{r}
positive.only <- function(x) { ifelse (x<0, 0, x) }
dem.ne.m <- calc(dem.ne.m, positive.only)
names(dem.ne.m)
names(dem.ne.m) <- "elev.m"
summary(dem.ne.m)
dem.ne.m.sp$elev.m <- pmax(0, dem.ne.m.sp$elev.m)
dem.ne.m.df$elev.m <- pmax(0, dem.ne.m.df$elev.m)
summary(dem.ne.m.df)
```

Save the objects (1) data frame, (2) projected, (3) geographic coords:

```{r}
save(dem.ne.m.df, dem.ne.m.sp, dem.ne.ll, dem.ne.ll.bbox, file="./dem_ne_4km.RData")
```

Save as a raster for use in QGIS:

```{r}
writeRaster(dem.ne.m, "./dem_ne_prov_4km.tif", format="GTiff",
            overwrite=TRUE)
writeRaster(dem.ne.m.bbox, "./dem_ne_bbox_4km.tif", format="GTiff",
            overwrite=TRUE)
```