Radar Aeroecology Workshop
Working with low tilt data for mapping animals near the ground
Working with low tilt data for mapping animals near the ground
July 20, 2023
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)
}