#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