1 Part 1

1.1 Setup

The following are libraries needed to run the code for this practical. They will be installed if necessary:

reqPackages <- c('geosphere','dplyr','elevatr','sf','fields','rgdal','foreign','terra','imputeTS','tzdb')
get.packages <- reqPackages[!(reqPackages %in% installed.packages()[,"Package"])]
if(length(get.packages)>0) install.packages(get.packages,dependencies = TRUE)

library(geosphere)
library(dplyr)
library(sf)
library(rgdal)
library(foreign)
library(terra)
library(imputeTS)
library(elevatr)
library(fields)

These are user defined variables of radar characteristics, file paths, and other constants. Please modify them as needed. The example radar is the KLWX NEXRAD station near Washington, DC USA. I have set the pulse volume dimensions to Super Resolution (0.5 degree by 250 m).

##User defined radar characteristics
site<-"KLWX"
X_origin <- -77.477778   #longitude of radar location
Y_origin <- 38.975278   #latitiude of radar location
zone <- 18 #UTM zone of radar
antht <- 118.54 #height of radar antenna above sea level in meters  = base ht (ASL) + tower ht

##define file paths
OUT <- "C:/Users/jbule/Dropbox/RAW_2022" #directory of files
DEM <- "lwxdem1.tif" #name of DEM raster file

setwd (OUT)
##User defined constants
degrees <- 0 #initial degree for creating sample volumes
step <- 0.5 #azimuth increment in degress (i.e., beamwidth)
depth <- 250 #depth of range gates (i.e., gatewidth)
minrange <- 2000 #min start range (m) for creating sample volumes
maxrange <- 5000 #max end range (m) for creating sample volumes
theta <- 0.5 #elevation angle of beam in degrees
halfbeam <- 0.008464847 #half power radar beam width in radians (e.g., 0.485 degree for NEXRAD)

The following are global constants that should not need to be modified.

##Constants
DEG2RAD <- 0.017453293 #convert degrees to radians
rad2deg <- 57.295779513 #convert radians to degrees
earthe  <- 6378137  #earth equatorial radius in meters

##compute equivalent earth radius at the height of the radar antenna
earth <- (earthe*(1 - (0.0033493*(sin(Y_origin*DEG2RAD))**2))) #earth radius to sea level at radar latititude
eartha <- (earth+antht) #Earth radius at antenna height
earthr <- 4*eartha/3 #Equivalent earth radius in a standard atmosphere

##prepare sequences for drawing spatial boundaries of sample volumes
azimuths<-seq(0, 359.5, by=step)
ids<-rep(1:(((maxrange-minrange)+depth)/depth*360/step), each=5)
ranges<-((seq(minrange,maxrange,by = depth))*cos(theta*DEG2RAD))
range<-rep(seq(minrange,maxrange,by = depth), each=720)

Here are some functions that will be used later:

#' Rounds to any decimal or whole value
#' 
#' @param x  vector of values
#' @param dec the number/decimal you want to convert to. 
mround <- function(x,dec){ dec*round(x/dec) }

#' Calculates the percentage of beam blockage across a radar range
#' 
#' This function requires input of a dataframe with three columns 'range','azimuth', and 'groundht'
#' @param data  data.frame with numberic columns for range, azimuth, and ground height (groundht)
#' @param anht the height of the radar antenna above sea level in m
#' @param theta  the tilt angle of the beam in degrees

beam_block<-function(data, theta, anht){
  
  data <- basegrid
  
  beam<-beamht_std(data$range, antht=antht, theta=theta, earthr=earthr)
  beam[,c('bbeam','mbeam','tbeam')] <- beam[,c('bbeam','mbeam','tbeam')] - data$groundht
  beam$block<-0
  beam$block[beam$tbeam<=0]<-1
  beam$block[beam$bbeam<0 & beam$mbeam >=0 ] <-mround(-1*beam$bbeam[beam$bbeam<0 & beam$mbeam >=0 ]/(2*beam$range[beam$bbeam<0 & beam$mbeam >=0 ]*halfbeam*cos(DEG2RAD*theta)), 0.05)
  beam$block[beam$bbeam <0 & beam$mbeam <0 ] <- mround(1-(beam$tbeam[beam$bbeam <0 & beam$mbeam <0 ]/(2*beam$range[ beam$bbeam <0 & beam$mbeam <0 ]*halfbeam*cos(DEG2RAD*theta))), 0.05)
  beam$azimuth <- data$azimuth
  beam <- beam[order(beam$azimuth, beam$range),]
  list1<-tapply(beam$block, beam$azimuth, list)
  
  here<-lapply(list1,  function(x) {
    listblock<-list()
    for(i in 1:length(x)){
      if(i==1){pblock<-0}
      pblock<-max(c(x[i],pblock))
      if(pblock<0){pblock<-0}
      listblock[[i]]<-pblock
    }
    listblock
  })
  
  
  beam$pblock<-unlist(here)
  beam
}

#' Calculates the top, middle, and bottom of a radar beam in a standard atmosphere
#'
#'  @param range a vector of ranges in m 
#'  @param earhtr the equivalent earth radius at the radar antenna in m
#'  @param anht the height of the radar antenna above sea level in m
#'  @param theta  the tilt angle of the beam in degrees

beamht_std<-function(range, earthr=earthr, antht=antht,theta=theta){
  
  mbeam <- sqrt(earthr^2+range^2+(2*range*earthr*sin(DEG2RAD*theta)))-(earthr)+antht #middle height of beam above sea level
  bbeam <- mbeam-range*halfbeam*(cos(DEG2RAD*0.5)) #bottom height of beam above sea level
  tbeam <- mbeam+(range*halfbeam*(cos(DEG2RAD*0.5)))#top height of beam above sea level
  data.frame(range=range,mbeam=mbeam, bbeam=bbeam, tbeam=tbeam, theta=theta)
  
}