horizon_4Oscar.R

library(sp)
library(raster)
library(rgdal)
library(maptools)
library(gstat)
library(lattice)
library(latticeExtra)
library(rasterVis)
library(parallel)
library(solaR)

## Configurate Projections
projUTM  <-  CRS("+proj=utm +zone=33 +ellps=WGS84")
projLongLat <- CRS(" +proj=longlat +ellps=WGS84")

# Load data 
load("elev_area.RData")
load("elev_Frame.RData")

## maximum sample distance
d <- 5000

## Sample distance defined by DEM resolution
resD <- max(res(elev))

## Separations
seps <- seq(resD, d, by=resD)

## Raster definition with UTMX, UTMY coordinates and elevation
xyelev <- stack(init(elev, v='x'),
                init(elev, v='y'),
                elev)
names(xyelev) <- c('x', 'y', 'elev')

## Coordinates
locs <- as.matrix(xyelev)
px <- locs[, 1]
py <- locs[, 2]
pz <- locs[, 3]

## Angle of sector sampling 1º
inc <- pi/180
alfa <- seq(-pi, pi - inc, inc)

## Compute the horizon angle for each direction (alfa) with parallel
## computing using mclapply.
old <- setwd(tempdir())
for (i in seq_along(alfa))
{
    ang <- alfa[i]
    cat("Angle: ", ang, "\n")
    ## Loop across separations
    lHor <- mclapply(seps, function(sep)
    {
        cat('Sep: ', sep, '\n')
        x1 <- px - sin(ang) * sep
        y1 <- py - cos(ang) * sep    
        p1 <- cbind(x1,y1)
        z1 <- extract(elevFrame, p1)
        ## Angle in degrees
        r2d(atan2(z1 - pz, sep))
    }, mc.cores = detectCores())
    ## The result is a list with an element for each separation. Each
    ## element is a vector of angles, with a cell for each point of
    ## the DEM (row of locs). The horizon angle at each location is
    ## the maximum angle, that must be positive. pmax computes the
    ## maximum value at each point.
    hor <- do.call(pmax, lHor)
    hor[hor < 0] <- 0
    ## Write the results to a raster file and remove the intermediate results
    r <- raster(elev)
    r <- setValues(r, hor)
    writeRaster(r, filename = paste0('horizon_', i), overwrite = TRUE)
    rm(r, lHor, hor)
    gc()
}

## The result is a set of files (RasterLayer), one per each direction angle.
hFiles <- paste0('horizon_', seq_along(alfa))
## Read them and build a RasterStack with a layer for each angle
horizon <- stack(hFiles)
writeRaster(horizon, 'horizonStack')
setwd(old)