R/horizon_1_1.R

In horizon: Horizon Search Algorithm

horizonSearch <- function(x, azimuth, maxDist=1000, degrees=FALSE, ll=TRUE, filename="", blockSize=NULL){
  if (degrees) azimuth <- azimuth*pi/180
  if (!nchar(filename)) filename <- rasterTmpFile()
  if (is.null(blockSize)) blockSize <- ceiling(nrow(x)/10)

  #Define size of buffer (from resolution in Y (in m) at UL/LR image)
  resY <- min(pointDistance(xyFromCell(x,cellFromRowCol(x,c(1,nrow(x)),c(1,ncol(x)))),
                            xyFromCell(x,cellFromRowCol(x,c(2,nrow(x)-1),c(1,ncol(x)))),
                            lonlat=ll))
  bufferSize <- ceiling(ceiling((maxDist/resY)*abs(cos(azimuth)))*1.05)*sign(cos(azimuth)) #(5% extra added for safety)
  
  #Normalized distance in X and Y
  dX <- sin(azimuth)/max(abs(c(sin(azimuth),cos(azimuth))))
  dY <- cos(azimuth)/max(abs(c(sin(azimuth),cos(azimuth))))

  #Initiate output raster
  outImg <- raster(x)
  outImg <- writeStart(outImg, filename, overwrite=TRUE, dataType="FLT4S")
  
  #Process each block of rows
  for(firstRow in seq(1,nrow(x),blockSize)){
    firstCell <- 1+(firstRow-1)*ncol(x)
    lastCell <- firstCell+blockSize*ncol(x)-1
    if (lastCell>ncell(x)) lastCell <- ncell(x)
    cells_row <- firstCell:lastCell
    elev_row <- getValues(x, firstRow, blockSize)

    firstRow_buffer <- max(1,min(firstRow,firstRow-bufferSize))
    lastRow_buffer <-min(max(firstRow+blockSize-1-bufferSize, firstRow+blockSize-1), nrow(x))
    firstCell_buffer <- 1+(firstRow_buffer-1)*ncol(x)
    lastCell_buffer <- lastRow_buffer*ncol(x)
    cells_buffer <- firstCell_buffer:lastCell_buffer
    elev_buffer <- getValues(x, firstRow_buffer, lastRow_buffer-firstRow_buffer+1)
    
    i=1
    cells_d <- cellFromRowCol(x,rowFromCell(x,cells_row)-round(i*dY),colFromCell(x,cells_row)+round(i*dX))
    valid <- validCell(x, cells_d)
    elev_d <- numeric(length=length(cells_d))+NA
    elev_d[valid] <- elev_buffer[cells_d[valid]-firstCell_buffer+1]
    
    dist_d <- pointDistance(xyFromCell(x,cells_row), xyFromCell(x,cells_d), lonlat=ll)
    meanDist <- mean(dist_d[valid],na.rm=TRUE)
    iHoriz <- pmax(0,atan2((elev_d - elev_row),dist_d))
    
    outHoriz_p <- numeric(length(cells_d))+NA
    valid <- valid & !is.na(iHoriz)
    outHoriz_p[valid] <- iHoriz[valid]

    while((meanDist < maxDist) & (sum(valid) > 0)){
      i=i+1
      cells_d <- cellFromRowCol(x,rowFromCell(x,cells_row)-round(i*dY),colFromCell(x,cells_row)+round(i*dX))
      valid <- validCell(x, cells_d)
      elev_d <- numeric(length=length(cells_d))+NA
      elev_d[valid] <- elev_buffer[cells_d[valid]-firstCell_buffer+1]
      
      dist_d <- pointDistance(xyFromCell(x,cells_row), xyFromCell(x,cells_d), lonlat=ll)
      meanDist <- mean(dist_d[valid],na.rm=TRUE)
      iHoriz <- pmax(0,atan2((elev_d - elev_row),dist_d))
      
      valid <- valid & !is.na(iHoriz) & (is.na(outHoriz_p) | (!is.na(iHoriz) & !is.na(outHoriz_p) & iHoriz > outHoriz_p))
      outHoriz_p[valid] <- iHoriz[valid]
    }
    if (degrees) outImg <- writeValues(outImg,outHoriz_p*180/pi, firstRow) else outImg <- writeValues(outImg,outHoriz_p, firstRow)
  }
  
  outImg <- writeStop(outImg)
  return(outImg)
}

svf <- function(x, nAngles=16, maxDist=1000, ll=TRUE, filename="", blockSize=NULL, verbose=TRUE){
  if (!nchar(filename)) filename <- rasterTmpFile()
  if (is.null(blockSize)) blockSize <- ceiling(nrow(x)/10) 
  angles <- seq(0,360-360/nAngles,360/nAngles)*pi/180
  
  #Derive terrain slope and azimuth
  sa <- terrain(x, opt=c("slope","aspect"),unit="radians")
  
  #Initiate output raster
  outImg <- raster(x)
  outImg <- writeStart(outImg, filename, overwrite=TRUE, dataType="FLT4S")
   
  resY <- min(pointDistance(xyFromCell(x,cellFromRowCol(x,c(1,nrow(x)),c(1,ncol(x)))),
                            xyFromCell(x,cellFromRowCol(x,c(2,nrow(x)-1),c(1,ncol(x)))),
                            lonlat=ll))
  
  if(verbose){
    nSteps <- ceiling(nrow(x)/blockSize)*nAngles
    step.count <- 0
    step.message <- seq(10,110,5)
    cat("0%...")
  }

  #Process each block of rows
  for(firstRow in seq(1,nrow(x),blockSize)){

    sl_as_row <- getValues(sa, firstRow, blockSize)
    svf_row <- numeric(nrow(sl_as_row))+0
  
    #Calculate horizon/sky view factor for each of nAngles
    for (azimuth in angles){
      if (verbose){
        step.count <- step.count+1
        step.frac <- 100*step.count/nSteps
        if(step.frac>step.message[1]){
          cat(step.message[1],"%...",sep="")
          step.message <- step.message[-1]
        }
      }

      #Define size of buffer (from resolution in Y (in m) at UL/LR image)
      bufferSize <- ceiling(ceiling((maxDist/resY)*abs(cos(azimuth)))*1.25)*sign(cos(azimuth)) #(25% extra added for safety)
    
      #Normalized distance in X and Y
      dX <- sin(azimuth)/max(abs(c(sin(azimuth),cos(azimuth))))
      dY <- cos(azimuth)/max(abs(c(sin(azimuth),cos(azimuth))))
    
      #Load elevation row + buffer
      firstCell <- 1+(firstRow-1)*ncol(x)
      lastCell <- firstCell+blockSize*ncol(x)-1
      if (lastCell>ncell(x)) lastCell <- ncell(x)
      cells_row <- firstCell:lastCell
      elev_row <- getValues(x, firstRow, blockSize)
      
      firstRow_buffer <- max(1,min(firstRow,firstRow-bufferSize))
      lastRow_buffer <-min(max(firstRow+blockSize-1-bufferSize, firstRow+blockSize-1), nrow(x))
      firstCell_buffer <- 1+(firstRow_buffer-1)*ncol(x)
      lastCell_buffer <- lastRow_buffer*ncol(x)
      cells_buffer <- firstCell_buffer:lastCell_buffer
      elev_buffer <- getValues(x, firstRow_buffer, lastRow_buffer-firstRow_buffer+1)
      
      i=1
      cells_d <- cellFromRowCol(x,rowFromCell(x,cells_row)-round(i*dY),colFromCell(x,cells_row)+round(i*dX))
      valid <- validCell(x, cells_d)
      elev_d <- numeric(length=length(cells_d))+NA
      elev_d[valid] <- elev_buffer[cells_d[valid]-firstCell_buffer+1]
      
      dist_d <- pointDistance(xyFromCell(x,cells_row), xyFromCell(x,cells_d), lonlat=ll)
      meanDist <- mean(dist_d[valid],na.rm=TRUE)
      iHoriz <- pmax(0,atan2((elev_d - elev_row),dist_d))
      
      outHoriz_p <- numeric(length(cells_d))+NA
      valid <- valid & !is.na(iHoriz)
      outHoriz_p[valid] <- iHoriz[valid]
      
      while((meanDist < maxDist) & sum(valid) > 0){
        i=i+1
        cells_d <- cellFromRowCol(x,rowFromCell(x,cells_row)-round(i*dY),colFromCell(x,cells_row)+round(i*dX))
        valid <- validCell(x, cells_d)
        elev_d <- numeric(length=length(cells_d))+NA
        elev_d[valid] <- elev_buffer[cells_d[valid]-firstCell_buffer+1]
        
        dist_d <- pointDistance(xyFromCell(x,cells_row), xyFromCell(x,cells_d), lonlat=ll)
        meanDist <- mean(dist_d[valid],na.rm=TRUE)
        iHoriz <- pmax(0,atan2((elev_d - elev_row),dist_d))
        
        valid <- valid & !is.na(iHoriz) & (is.na(outHoriz_p) | (!is.na(iHoriz) & !is.na(outHoriz_p) & iHoriz > outHoriz_p))
        outHoriz_p[valid] <- iHoriz[valid]
      }
      
      #Correct horizon angle in selected azimuth based on 8-neighbour slope
      horMin <- -atan(tan(sl_as_row[,1])*cos(azimuth-sl_as_row[,2]))
      outHoriz_p <- pmax(horMin, outHoriz_p, na.rm=FALSE)
      
      #Calculate SVF in angle direction from horizon + terrain slope/aspect
      svf_angle <- ((cos(sl_as_row[,1])*(cos(outHoriz_p)^2))+
                      sin(sl_as_row[,1])*cos(azimuth-sl_as_row[,2])*((pi/2)-outHoriz_p-(sin(outHoriz_p)*cos(outHoriz_p))))
      svf_row <- svf_row+svf_angle/nAngles

    }
    outImg <- writeValues(outImg,svf_row,firstRow)
  }
  outImg <- writeStop(outImg)
  return(outImg) 
}