Skip to content

Instantly share code, notes, and snippets.

@oscarperpinan
Created June 20, 2016 16:50
Show Gist options
  • Save oscarperpinan/f9933e43447e2de233f7ab95e7ae1c22 to your computer and use it in GitHub Desktop.
Save oscarperpinan/f9933e43447e2de233f7ab95e7ae1c22 to your computer and use it in GitHub Desktop.
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment