Shaded relief by clustering aspect

#Shaded relief by clustering aspect
#Fabio Veronesi

#Loading packages
library(raster)
library(cluster)
library(RSAGA)

#Working directory.
setwd("...")


###################
###  FUNCTIONS  ###
###################
#Clustering
Clustering <- function(grid){
values<-getValues(grid)
i <- which(! is.na(values))
clus<-clara(values[i],4)
clusters<-raster(aspect)
clusters[i]<-clus$clustering
writeRaster(clusters,filename="Clusters.asc",format="ascii")}

#Majority filter
Majority_Filter <- function(grid,radius){
rsaga.esri.to.sgrd(in.grids=grid@file@name,"Clusters.sgrd")
rsaga.geoprocessor(lib="grid_filter",module=6,list(INPUT="Clusters.sgrd",RESULT="Clust_Maj.sgrd",MODE=0,RADIUS=radius))
rsaga.sgrd.to.esri("Clust_Maj.sgrd","Clust_Maj.asc")
raster("Clust_Maj.asc")}

#Mean filter
Mean_Filter <- function(grid,radius){
#rsaga.esri.to.sgrd(in.grids=grid@file@name,"Clust_Maj.sgrd")
rsaga.geoprocessor(lib="geostatistics_grid",module=1,list(GRID="Clust_Maj.sgrd",MEAN="Clust_Mean.sgrd",DIFF="tmp.sgrd",STDDEV="tmp.sgrd",RANGE="tmp.sgrd",MIN="tmp.sgrd",MAX="tmp.sgrd",DEVMEAN="tmp.sgrd",PERCENT="tmp.sgrd",RADIUS=radius))
rsaga.sgrd.to.esri("Clust_Mean.sgrd","Clust_Mean.asc")
raster("Clust_Mean.asc")}

#Sine-Wave
sine.wave<-function(x){((MAX.LIGHT-MIN.LIGHT)/2)*sin(3.141593*x-1.570796)+(MAX.LIGHT-((MAX.LIGHT-MIN.LIGHT)/2))}


#####################
###  USER INPUTS  ###
#####################
MAX.LIGHT = 315
MIN.LIGHT = 280

MAJ.F.RADIUS = 25
MEA.F.RADIUS = 25

ZENITH.OPTION = 1
#[0] = Constant zenith (please set ZENITH.VALUE)
#[1] = Zenith changes with Elevation
#[2] = Zenith changes with Slope

ZENITH.VALUE = 45
#Used only if ZENITH.OPTION is equal to 0



###############
###  START  ###
###############

#Loading DTM
dtm<-raster("dtm_10m.asc")
proj4string(dtm)=CRS("+init=epsg:2056")



#Derivatives
slope<-terrain(dtm,opt='slope',unit='radians',neighbors=8)
aspect<-terrain(dtm,opt='aspect',unit='degrees',neighbors=8)



#Clustering
cluster_raster <- Clustering(aspect)



#Majority Filter
majority <- Majority_Filter(cluster_raster,MAJ.F.RADIUS)


#Mean Filter
mean <- Mean_Filter(majority,MEA.F.RADIUS)



#Azimuth
azimuth<-raster(dtm)
names(azimuth)<-"azimuth"
azimuth[]<-sine.wave(mean[])*(pi/180)



#Zenith
if(ZENITH.OPTION==0){
#Constant Zenith
zenith<-raster(dtm)
names(zenith)<-"zenith"
zenith[]<- rep(ZENITH.VALUE,length(dtm[]))
}

if(ZENITH.OPTION==1){
#Zenith Elevation
zenith<-raster(dtm)
names(zenith)<-"zenith"
slo<-(dtm[]-min(dtm[],na.rm=T))/(max(dtm[],na.rm=T)-min(dtm[],na.rm=T))
zenith[]<-(70*slo)*(pi/180)
}

if(ZENITH.OPTION==2){
#Zenith Slope
zenith<-raster(slope)
names(zenith)<-"zenith"
slo<-(slope[]-min(slope[],na.rm=T))/(max(slope[],na.rm=T)-min(slope[],na.rm=T))
zenith[]<-(70*slo)
}




#Hillshading
aspect<-terrain(dtm,opt='aspect',unit='radians',neighbors=8)
staked<-stack(slope,aspect,azimuth,zenith)
hill.shd<-raster(aspect)
Shading<-((cos(getValues(staked$zenith))*cos(getValues(staked$slope)))+(sin(getValues(staked$zenith))*sin(getValues(staked$slope))*cos(getValues(staked$azimuth)-getValues(staked$aspect))))
Shading[Shading<0]<-0
hill.shd[]<-255*Shading



#View Results
image(hill.shd,col=gray.colors(255))


#Save in ASC Format