Calculate slope

slope=terrain(reg1,opt="slope",unit="degrees")
plot(slope)

Calculate aspect

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.

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

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()

TRI (Terrain Ruggedness Index)

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

tri=terrain(reg1,opt="TRI")
plot(tri)

Roughness

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

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.

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:

flowdir=terrain(reg1,opt="flowdir")

plot(flowdir)

Much more powerful hydrologic modeling in GRASS GIS

Gridded Population of the World

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

The steps to aquire the full dataset are as follows:

• Log into SEDAC with an Earth Data Account 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.

pop_global=raster(file.path(datadir,"gpw-v4-population-density-2015/gpw-v4-population-density_2015.tif"))

data(bangladesh_population)

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

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

## make a virtual copy with a shorter name for convenience

pop=bangladesh_population

Explore population data

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.

pop

## class       : RasterLayer 
## dimensions  : 707, 560, 395920  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 88.00833, 92.675, 20.74167, 26.63333  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : gpw.v4.population.density_2015 
## values      : 19.83075, 154258.4  (min, max)

dem

## class       : RasterLayer 
## dimensions  : 3000, 1200, 3600000  (nrow, ncol, ncell)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : 89.99958, 90.99958, 21.49958, 23.99958  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/adamw/Downloads/data/dem_bgd.tif 
## names       : dem_bgd 
## values      : -27, 37  (min, max)

res(pop)

## [1] 0.008333333 0.008333333

res(dem)

## [1] 0.0008333333 0.0008333333

origin(pop)

## [1] -6.536993e-13  2.629008e-13

origin(dem)

## [1] -0.000416061 -0.000416207

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

## [1] 0.7828593

cellStats(raster::area(dem),"mean")

## [1] 0.007886292

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

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()

floodpop2=flood2_area*pop_fine
cellStats(floodpop2,sum)

## [1] 29796929

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()

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

res(pop)/res(dem)

## [1] 10 10

dem_coarse=dem%>%
  aggregate(fact=10,fun=min,expand=T)%>%
  resample(pop,method="bilinear")

For the coarse resolution data

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)

## [1] 40077457