06_CreatingWorkflows.Rmd

---
title: "Creating Workflows"

---

```{r, echo=FALSE, message=FALSE, results='hide', purl=FALSE}
source("knitr_header.R")
```

[<i class="fa fa-file-code-o fa-3x" aria-hidden="true"></i> The R Script associated with this page is available here](`r output`).  Download this file and open it (or copy-paste into a new script) with RStudio so you can follow along.  


## Libraries

```{r message=F,warning=FALSE}
library(knitr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(raster)
library(rasterVis)
library(scales)
library(rgeos)

# load data for this course
# devtools::install_github("adammwilson/DataScienceData")
library(DataScienceData)

```

## Today's question

### How will future (projected) sea level rise affect Bangladesh?

1. How much area is likely to be flooded by rising sea level?
2. How many people are likely to be displaced?
3. Will sea level rise affect any major population centers?

## Bangladesh

```{r}
getData("ISO3")%>%
  as.data.frame%>%
  filter(NAME=="Bangladesh")
```




### Download Bangladesh Border

Often good idea to keep data in separate folder.  You will need to edit this for your machine!
```{r}
datadir="~/Downloads/data"
if(!file.exists(datadir)) dir.create(datadir, recursive=T)
```
Download country border.
```{r, eval=F}
bgd=getData('GADM', country='BGD', level=0,path = datadir)
```

Or load it from the data package.
```{r}
data(bangladesh)
bgd=bangladesh
```

```{r}
bgd%>%
  gSimplify(0.01)%>%
  plot()
```

## Topography

SRTM Elevation data with `getData()` as 5deg tiles.  If you have trouble downloading using `getData()`, skip to the `data(bangladesh_dem)` line below

```{r, eval=F}
bgdc=gCentroid(bgd)%>%coordinates()

dem1=getData("SRTM",lat=bgdc[2],lon=bgdc[1],path=datadir)
```


### Mosaicing/Merging rasters

Download the remaining necessary tiles

```{r, eval=F}
dem2=getData("SRTM",lat=23.7,lon=85,path=datadir)
```

Use `merge()` to join two aligned rasters (origin, resolution, and projection).  Or `mosaic()` combines with a function.

```{r, eval=F}
dem=merge(dem1,dem2)
```

Or, load it from the data package.
```{r}
data(bangladesh_dem)
dem=bangladesh_dem # rename for convenience
```

```{r}
plot(dem)
bgd%>%
  gSimplify(0.01)%>%
  plot(add=T)
```

## Saving/exporting rasters

Beware of massive temporary files!

```{r}
inMemory(dem)
dem@file@name
file.size(sub("grd","gri",dem@file@name))*1e-6
showTmpFiles()
```



```{r}
rasterOptions()
```
Set with `rasterOptions(tmpdir = "/tmp")`



Saving raster to file: _two options_

Save while creating
```{r,eval=F}
dem=merge(dem1,dem2,filename=file.path(datadir,"dem.tif"),overwrite=T)
```

Or after
```{r,eval=F}
writeRaster(dem, filename = file.path(datadir,"dem.tif"))
```



### WriteRaster formats

Filetype  Long name	                      Default extension	  Multiband support
---       ---                             ---                 ---                                                     
raster	  'Native' raster package format	.grd	              Yes
ascii	    ESRI Ascii	                    .asc                No
SAGA	    SAGA GIS	                      .sdat	              No
IDRISI	  IDRISI	                        .rst	              No
CDF	      netCDF (requires `ncdf`)	      .nc	                Yes
GTiff	    GeoTiff (requires rgdal)	      .tif	              Yes
ENVI	    ENVI .hdr Labelled	            .envi	              Yes
EHdr	    ESRI .hdr Labelled	            .bil	              Yes
HFA	      Erdas Imagine Images (.img)   	.img	              Yes

`rgdal` package does even more...

### Crop to Coastal area of Bangladesh

```{r}
# crop to a lat-lon box
dem=crop(dem,extent(90,91,21.5,24),filename=file.path(datadir,"dem_bgd.tif"),overwrite=T)

plot(dem)
bgd%>%
  gSimplify(0.01)%>%
  plot(add=T)
```

# Use ggplot
```{r,warning=F}
gplot(dem,max=1e5)+
  geom_tile(aes(fill=value))+
  scale_fill_gradientn(
    colours=c("red","yellow","grey30","grey20","grey10"),
    trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e3)))+
  coord_equal(ylim=c(21.5,24),xlim=c(90,91))+
  geom_path(data=fortify(bgd),
            aes(x=long,y=lat,group=group),size=.5)
```



# Terrain analysis (an aside)

## Terrain analysis options

`terrain()` options:

* slope
* aspect
* TPI (Topographic Position Index)
* TRI (Terrain Ruggedness Index)
* roughness
* flowdir



Use an even smaller region:
```{r}
reg1=crop(dem,extent(90.6,90.7,23.25,23.4))
plot(reg1)
```

The terrain indices are according to Wilson et al. (2007), as in [gdaldem](http://www.gdal.org/gdaldem.html).



### Calculate slope

```{r}
slope=terrain(reg1,opt="slope",unit="degrees")
plot(slope)
```



### Calculate aspect

```{r}
aspect=terrain(reg1,opt="aspect",unit="degrees")
plot(aspect)
```



### TPI (Topographic Position Index)

Difference between the value of a cell and the mean value of its 8 surrounding cells.

```{r}
tpi=terrain(reg1,opt="TPI")

gplot(tpi,max=1e6)+geom_tile(aes(fill=value))+
  scale_fill_gradient2(low="blue",high="red",midpoint=0)+
  coord_equal()
```
Negative values indicate valleys, near zero flat or mid-slope, and positive ridge and hill tops



<div class="well">
## Your turn

* Identify all the pixels with a TPI less than -5 or greater than 5.
* Use `plot()` to:
    * plot elevation for this region
    * overlay the valley pixels in blue
    * overlay the ridge pixels in red

Hint: use `transparent` to plot a transparent pixel and `add=T` to add a layer to an existing plot. 

<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo1">Show Solution</button>
<div id="demo1" class="collapse">

```{r,purl=F}
plot(reg1)
plot(tpi>5,col=c("transparent","red"),add=T,legend=F)
plot(tpi<(-5),col=c("transparent","blue"),add=T,legend=F)

#OR (ggplot solution, sort of)
rcl=matrix(c(-Inf,-5,1,
           -5,5,2,
           5,Inf,3),byrow=T,nrow=3)
regclass=reclassify(tpi,rcl)
gplot(regclass,max=1e6)+geom_tile(aes(fill=value))+
  scale_fill_gradient2(low="blue",high="red",midpoint=2)+
  coord_equal()
```
</div>
</div>


### TRI (Terrain Ruggedness Index)

Mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells.

```{r}
tri=terrain(reg1,opt="TRI")
plot(tri)
```



### Roughness 

Difference between the maximum and the minimum value of a cell and its 8 surrounding cells.

```{r}
rough=terrain(reg1,opt="roughness")
plot(rough)
```




### Hillshade (pretty...)

Compute from slope and aspect (in radians). Often used as a backdrop for another semi-transparent layer.

```{r}
hs=hillShade(slope*pi/180,aspect*pi/180)

plot(hs, col=grey(0:100/100), legend=FALSE)
plot(reg1, col=terrain.colors(25, alpha=0.5), add=TRUE)
```



### Flow Direction

_Flow direction_ (of water), i.e. the direction of the greatest drop in elevation (or the smallest rise if all neighbors are higher). 

Encoded as powers of 2 (0 to 7). The cell to the right of the focal cell 'x' is 1, the one below that is 2, and so on:

32	64	    128
--- ---     ---     
16	**x**	  1
8   4       2



```{r}
flowdir=terrain(reg1,opt="flowdir")

plot(flowdir)
```
Much more powerful hydrologic modeling in [GRASS GIS](https://grass.osgeo.org) 

# Sea Level Rise



## Global SLR Scenarios

```{r,results="markdown"}
slr=data.frame(year=2100,
               scenario=c("RCP2.6","RCP4.5","RCP6.0","RCP8.5"),
               low=c(0.26,0.32,0.33,0.53),
               high=c(0.54,0.62,0.62,0.97))
slr
```

[IPCC AR5 WG1 Section 13-4](https://www.ipcc.ch/pdf/assessment-report/ar5/wg1/drafts/fgd/WGIAR5_WGI-12Doc2b_FinalDraft_Chapter13.pdf)

## Storm Surges

Range from 2.5-10m in Bangladesh since 1960 [Karim & Mimura, 2008](http://www.sciencedirect.com/science/article/pii/S0959378008000447).  

```{r}
ss=c(2.5,10)
```

## Raster area

1st Question: How much area is likely to be flooded by rising sea levels? 

WGS84 data is unprojected, must account for cell area (in km^2)...
```{r}
area=raster::area(dem)
plot(area)
```


<div class="well">
## Your Turn

1. How much area is likely to be flooded by rising sea levels for two scenarios:
   * 0.26m SLR and 2.5m surge (`r .26+2.5` total) - call this `flood1`
   * 0.97 SLR and 10m surge (`r 0.97+10` total) - call this `flood2`
   
Steps:

* Identify which pixels are below thresholds
* Multiply by cell area
* Use `cellStats()` to calculate potentially flooded areas.

<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo2">Show Solution</button>
<div id="demo2" class="collapse">
## Identify pixels below thresholds

```{r, purl=F}
flood1=dem<=2.76
flood2=dem<=10.97

plot(flood2,col=c("transparent","darkred"))
plot(flood1,col=c("transparent","red"),add=T)
```


## Multiply by area and sum

```{r, purl=F}
flood1_area=flood1*area
flood2_area=flood2*area

cellStats(flood1_area,sum)
cellStats(flood2_area,sum)
```

</div>
</div>

## Reclassification

Another useful function for raster processing is `reclass()`.

```{r}
rcl=matrix(c(-Inf,2.76,1,
           2.76,10.97,2,
           10.97,Inf,3),byrow=T,ncol=3)
rcl
regclass=reclassify(dem,rcl)

gplot(regclass,max=1e5)+
  geom_tile(aes(fill=as.factor(value)))+
  scale_fill_manual(values=c("red","orange","blue"),
                    name="Flood Class")+
  coord_equal()
```


Or, do reclassification 'on the fly in the plotting function

```{r}
gplot(dem,max=1e5)+
  geom_tile(aes(fill=cut(value,c(-Inf,2.76,10.97,Inf))))+
  scale_fill_manual(values=c("red","orange","blue"),
                    name="Flood Class")+
  coord_equal()
```



## Socioeconomic Data

Socioeconomic Data and Applications Center (SEDAC)
[http://sedac.ciesin.columbia.edu](http://sedac.ciesin.columbia.edu)
<img src="06_assets/sedac.png" alt="alt text" width="70%">

* Population
* Pollution
* Energy
* Agriculture
* Roads



### Gridded Population of the World

Data _not_ available for direct download (e.g. `download.file()`) and are only available globally.  

<img src="06_assets/sedacData.png" alt="alt text" width="80%">

The steps to aquire the full dataset are as follows:

* Log into SEDAC with an Earth Data Account
[http://sedac.ciesin.columbia.edu](http://sedac.ciesin.columbia.edu)
* Download Population Density Grid for 2015
* Crop and mask to the country boundary for Bangladesh

The masked data are available in the DataScienceData package in the `bangladesh_pop` dataset.

### Load population data

Use `raster()` to load a raster from disk.

```{r, eval=F}
pop_global=raster(file.path(datadir,"gpw-v4-population-density-2015/gpw-v4-population-density_2015.tif"))
```

```{r}
data(bangladesh_population)
```

If the data package isn't working, download directly from github.

```{r, eval=F}
tf=tempfile()
download.file("https://github.com/adammwilson/DataScienceData/raw/master/data/bangladesh_population.rda",destfile = tf)
load(tf)
```

```{r}
## make a virtual copy with a shorter name for convenience

pop=bangladesh_population
```

### Explore population data

```{r}
gplot(pop,max=1e5)+geom_tile(aes(fill=value))+
  scale_fill_gradientn(colours=c("grey90","grey60","darkblue","blue","red"),
                       trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e5)))+
  coord_equal()
```



### Resample to DEM

Compare the resolution and origin of `pop` and `dem`.

```{r}
pop
dem

res(pop)
res(dem)

origin(pop)
origin(dem)

# Look at average cell area in km^2 
cellStats(raster::area(pop),"mean")
cellStats(raster::area(dem),"mean")
```

So to work with these rasters (population and elevation), it is easiest to "adjust" them to have the same resolution.  But there is no good way to do this.  Do you aggregate the finer raster or resample the coarser one?

Assume equal density within each grid cell and resample
```{r, warning=F}
pop_fine=pop%>%
  resample(dem,method="bilinear")

gplot(pop_fine,max=1e5)+geom_tile(aes(fill=value))+
  scale_fill_gradientn(
    colours=c("grey90","grey60","darkblue","blue","red"),
    trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e5)))+
  coord_equal()

```


<div class="well">
## Your turn

How many people are likely to be displaced?

Steps:

* Multiply flooded area (`flood2`) **x** population density **x** area
* Summarize with `cellStats()`
* Plot a map of the number of people potentially affected by `flood2`

<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo3">Show Solution</button>
<div id="demo3" class="collapse">

For the fine resolution population data
```{r, purl=F}
floodpop2=flood2_area*pop_fine
cellStats(floodpop2,sum)
```


Number of potentially affected people across the region.

```{r, purl=F}
gplot(floodpop2,max=1e6)+geom_tile(aes(fill=value))+
  scale_fill_gradientn(
    colours=c("grey90","grey60","darkblue","blue","red"),
    trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e4)))+
  coord_equal()
```

</div>
</div>

Or resample elevation to resolution of population:
1. First aggregate to approximate spatial resolution
2. Resample to align grids perfectly

```{r}
res(pop)/res(dem)
dem_coarse=dem%>%
  aggregate(fact=10,fun=min,expand=T)%>%
  resample(pop,method="bilinear")
```

For the coarse resolution data
```{r, purl=F}
flood_coarse=dem_coarse<=10.97
dem_coarse_area=raster::area(dem_coarse)
flood_coarse_area=flood_coarse*dem_coarse_area
floodpop_coarse=flood_coarse_area*pop
cellStats(floodpop_coarse,sum)
```


## Raster Distances

`distance()` calculates distances for all cells that are NA to the nearest cell that is not NA.

```{r}
popcenter=pop>5000
popcenter=mask(popcenter,popcenter,maskvalue=0)
plot(popcenter,col="red",legend=F)
```


In meters if the RasterLayer is not projected (`+proj=longlat`) and in map units (typically also meters) when it is projected.

```{r, warning=F}
popcenterdist=distance(popcenter)
plot(popcenterdist)
```


<div class="well">
## Your Turn

Will sea level rise affect any major population centers?

Steps:

* Resample `popcenter` to resolution of `dem` using `method=ngb`
* Identify `popcenter` areas that flood according to `flood2`.


<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo4">Show Solution</button>
<div id="demo4" class="collapse">

Will sea level rise affect any major population centers?

```{r,purl=F}
popcenter2=raster::resample(popcenter,dem,method="ngb")

floodpop2= flood2==1 & popcenter2
floodpop2=mask(floodpop2,floodpop2,maskval=0)

plot(flood2);plot(floodpop2,add=T,col="red",legend=F);
bgd%>%
  gSimplify(0.01)%>%
  plot(add=T)
```

</div>
</div>

## Vectorize raster

```{r,warning=F, message=F}
vpop=rasterToPolygons(popcenter, dissolve=TRUE)

gplot(dem,max=1e5)+geom_tile(aes(fill=value))+
  scale_fill_gradientn(
    colours=c("red","yellow","grey30","grey20","grey10"),
    trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e3)))+
  coord_cartesian(ylim=c(21.5,24),xlim=c(90,91))+
  geom_path(data=fortify(bgd),aes(x=long,y=lat,group=group),size=.5)+
  geom_path(data=fortify(vpop),aes(x=long,y=lat,group=group),size=1,col="green")

```
Warning: very slow on large rasters...

## 3D Visualization
Uses `rgl` library.  

```{r, eval=F}
plot3D(dem)
decorate3d()
```

<img src="06_assets/plot3d.png" alt="alt text" width="70%">

50 different styles illustrated [here](https://cran.r-project.org/web/packages/plot3D/vignettes/volcano.pdf).



Overlay population with `drape`

```{r, eval=F}
plot3D(dem,drape=pop, zfac=0.2)
decorate3d()
```

## Raster overview

* Perform many GIS operations
* Convenient processing and workflows
* Some functions (e.g. `distance()`) can be slow!