find_ridges.R

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