Skip to content

Instantly share code, notes, and snippets.

@rmarrotte
Last active July 29, 2016 15:57
Show Gist options
  • Save rmarrotte/209f2e0f98933ebcee67e5f277b388e2 to your computer and use it in GitHub Desktop.
Save rmarrotte/209f2e0f98933ebcee67e5f277b388e2 to your computer and use it in GitHub Desktop.
# Clear workspace
rm(list=ls())
# Load libraries
require(rasterVis)
require(dismo)
require(GISTools)
require(spdep)
require(RColorBrewer)
# Load my Digital Elevation Model
dem <- getData('alt', country='CHE')
# Plot it
x11(12,12)
levelplot(dem,col.regions=terrain.colors(20),
margin = FALSE)
# Focal stats using a 3 by 3 moving window
mean_filter <- focal(dem,w=matrix(1/(3^2), nc=3, nr=3))
# Stack them so we can plot them together
x11(24,12)
levelplot(stack(dem,mean_filter),col.regions=terrain.colors(20),
layout=c(2, 1),
names.attr = c("Digital Elevation Model", "Average filter"),
margin = FALSE,
colorkey = list(space = "bottom"))
# Topographic position index
TPI <- dem - mean_filter
x11(12,12)
levelplot(TPI, main = "Topographic Position Index",
par.settings=rasterTheme(region=brewer.pal(11,'RdYlBu')),
margin = FALSE)
# Some stats
quantile(TPI[],na.rm=T)
hist(TPI[])
# Get the standard deviation of the TPI
SD <- sd(TPI[],na.rm=T)
# Make landform classes
#Morphologic class De Reu et al. 2013; Weiss (2001)
landform <- reclassify(TPI, matrix(c(-Inf, -SD, 1,
-SD, -SD/2, 2,
-SD/2, 0, 3,
0, SD/2, 4,
SD/2, SD, 5,
SD, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform <- as.factor(landform)
rat <- levels(landform)[[1]]
rat[["landcover"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform) <- rat
# Plot the classification
x11(12,12)
levelplot(landform, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat$landcover,
main = "Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat[["landcover"]])))
# Extract ridges from the classification
ridges <- landform
ridges[] <- ifelse(ridges[] != 6, NA, 6)
# Reset the levels
ridges <- ratify(ridges)
# Adjust extent and get GM terrain map
ridgesGM <- projectRaster(from = ridges, crs = CRS("+init=epsg:3857"),method = "ngb") # takes awhile > 5 mins
mygmap <- gmap(x = ridgesGM, type = "terrain", zoom = 9)
ridgesGM_c <- crop(ridgesGM, mygmap)
# Plot it
x11(16,12)
plot(mygmap)
plot(ridgesGM_c, col = adjustcolor("red", alpha.f = 0.4) , add = T, legend = F)
rect(xleft = xmax(ridgesGM_c)-5000, ybottom = ymax(ridgesGM_c)-30000,
xright = xmax(ridgesGM_c)-30000, ytop = ymax(ridgesGM_c)-1000,
col= adjustcolor(col = "white", alpha.f = 0.7),
border = NA)
map.scale(xc = xmax(ridgesGM_c)-18000, yc = ymax(ridgesGM_c)-22000, len = 10000, ndivs = 1, subdiv = 10, units = "Kilometers")
north.arrow(xb= xmax(ridgesGM_c)-18000, yb= ymax(ridgesGM_c)-13000, len = 2000, lwd = 2)
legend(x = xmin(ridgesGM_c)+5000, y = ymax(ridgesGM_c)-1000,legend = "Ridges",
fill = adjustcolor("red", alpha.f = 0.4), inset = c(0,0),
bg = adjustcolor("white", alpha.f = 0.7), box.lty=0)
# Convert the ridges to points
ridgesGM_pnt <- rasterToPoints(ridgesGM_c, spatial = T)
# Make neighbours 2 km
edges <- dnearneigh(coordinates(ridgesGM_pnt), d1 = 0, d2 = 2000, longlat = F, row.names = row.names(ridgesGM_pnt))
# Add in how many links per node to the sp obj
ridgesGM_pnt$Nber_con <- as.numeric(lapply(edges, function(X){length(X)}))
table(ridgesGM_pnt$Nber_con) # check it
# Create a function to generate a continuous color palette
rbPal <- adjustcolor(rev(brewer.pal(4,"RdBu")),alpha.f = 0.7)
# Add it to the data
ridgesGM_pnt$Col <- rbPal[as.numeric(cut(ridgesGM_pnt$Nber_con,breaks = 4))]
# Plot the graph
x11(16,12)
plot(mygmap)
points(ridgesGM_pnt, col = ridgesGM_pnt$Col, pch = 19)
rect(xleft = xmax(ridgesGM_c)-5000, ybottom = ymax(ridgesGM_c)-30000,
xright = xmax(ridgesGM_c)-30000, ytop = ymax(ridgesGM_c)-1000,
col= adjustcolor(col = "white", alpha.f = 0.7),
border = NA)
map.scale(xc = xmax(ridgesGM_c)-18000, yc = ymax(ridgesGM_c)-22000, len = 10000, ndivs = 1, subdiv = 10, units = "Kilometers")
north.arrow(xb= xmax(ridgesGM_c)-18000, yb= ymax(ridgesGM_c)-13000, len = 2000, lwd = 2)
legend(x = xmin(ridgesGM_c)+5000, y = ymax(ridgesGM_c)-1000,
title = "# links",
legend = c("[0-2]","[3-4]","[5-6]","[7-8]"),
pch = 19, col = rbPal, inset = c(0,0),
bg = adjustcolor("white", alpha.f = 0.9), box.lty=0)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment