R/RNomadsTools.R

Defines functions SubsetNOMADS PlotWindProfile LinkExtractor MagnitudeAzimuth ModelGrid BuildProfile GetClosestForecasts

Documented in BuildProfile GetClosestForecasts LinkExtractor MagnitudeAzimuth ModelGrid PlotWindProfile SubsetNOMADS

#Functions for doing specific useful tasks with rNOMADS

GetClosestForecasts <- function(abbrev, forecast.date, model.date = "latest", depth = NULL, verbose = TRUE) {
 #Figure out the closest forecast to a given date, returns both the closest forecast behind and the closest forecast ahead, as well as how far behind and how far ahead
   #INPUTS
   #    ABBREV - Which model product you want
   #    FORECAST.DATE - What date you want a forecast for, as a date/time object, in GMT
   #    MODEL.DATE - Which model run to use, in YYYYMMDDHH, where HH is 00, 06, 12, 18.  Defaults to "latest", which means get the latest model available.
   #    DEPTH - How many model URLS to exame.  This is only taken into account when model.date!="latest".
   #       The option allows users to prevent the code from scanning every single URL.
   #        If MODEL.DATE is not "latest", rNOMADS will scan the entire directory, and this will take time. 
   #    VERBOSE - Give a blow-by-blow account of progress.  Defaults to TRUE.
   #OUTPUTS
   #    MODEL.URL - Which model URL to use
   #    MODEL.RUN.DATE - When the model was run
   #    BACK.FORECAST - Nearest forecast string behind forecast date
   #    FORE.FORECAST - Nearest forecast past forecast date
   #    BACK.HR - How many hours before the requested time the back forecast is
   #    FORE.HR - How many hours after the requested time the forward forecast is

   if(verbose){
       print("Finding model run dates...")
   }
   if(model.date == "latest")  {
       urls.out <- CrawlModels(abbrev = abbrev, depth = 2, verbose = verbose)
       model.parameters <- ParseModelPage(urls.out[1])
       if(length(model.parameters$pred) == 0) { #No data in latest model, try next latest
           model.parameters <- ParseModelPage(urls.out[2]) 
           if(length(model.parameters$pred) == 0) { #Nothing here either, maybe the server is down?
               stop("No data found for most recent model run.  Perhaps the NOMADS server is down?")
           } else {
               url.to.use <- urls.out[2]
           }
       } else {
           url.to.use <- urls.out[1]
       }
    } else {
         urls.out <- CrawlModels(abbrev = abbrev, depth = depth, verbose = verbose)
         model.run.dates <- unlist(stringr::str_extract_all(urls.out, "\\d{10}")) 
         if(model.date %in% model.run.dates) {
              url.to.use <- urls.out[which(model.date == model.run.dates)]
              model.parameters <- ParseModelPage(url.to.use)
         } else {
             stop(paste0("The model run date ", model.date, " does not appear in the list of model runs on the NOMADS server."))
         }
    }

    if(verbose) {
        print("Determining closest forecasts...")
    }
   
    #Get model run date, convert to POSIX date 
    run.date <- stringr::str_match_all(url.to.use, "\\d{10}")[[1]][1,1]
    d.vec <- strsplit(run.date, split = "")[[1]]
    nice.run.date <- strftime(paste0(paste(d.vec[1:4], collapse = ""),
        "-", paste(d.vec[5:6], collapse = ""), "-", paste(d.vec[7:8], collapse = ""),
        " ", paste(d.vec[9:10], collapse = ""), ":00:00", sep = "")) 

   #Figure out time difference between forecast date and model date 
   hr.shift <- as.numeric(difftime(forecast.date, as.POSIXlt(nice.run.date, tz = "GMT")), units = "hours")

   #Figure out the forward (in the future) forecast and back (in the past) forecast using hr.shift
   pred.hrs <- as.numeric(unlist(stringr::str_match_all(model.parameters$pred, "\\d{2,3}$")))
   hr.diff <- pred.hrs - hr.shift

   #Get forecasts
   back.forecast <- model.parameters$pred[which(max(hr.diff[which(hr.diff <=0)]) == hr.diff)]
   fore.forecast <- model.parameters$pred[which(min(hr.diff[which(hr.diff > 0)]) == hr.diff)]
   back.hr <- hr.diff[which(max(hr.diff[which(hr.diff <=0)]) == hr.diff)]
   fore.hr <- hr.diff[which(min(hr.diff[which(hr.diff > 0)]) == hr.diff)] 

   if(verbose) {
       print("Finished determining forecast times.")
   }
   return(list(model.url = url.to.use, model.run.date = nice.run.date, back.forecast = back.forecast, fore.forecast = fore.forecast, back.hr = back.hr, fore.hr = fore.hr))
}

BuildProfile <- function(model.data, lon, lat, spatial.average = FALSE, points = 4) {
    #This function builds an atmospheric profile, performing spatial interpolation if requested
    #INPUTS
    #    MODEL.DATA - Data structure returned from ReadGrib or DODSGrab
    #    LON - Longitude of points of interest
    #    LAT - Latitude of points of interest
    #    SPATIAL.AVERAGE - Boolean determining whether to get nearest node value (FALSE) or interpolate using b-splines (TRUE)
    #    POINTS - If SPATIAL.AVERAGE = TRUE, use this many points proximal to the point of interest to interpolate
    #OUTPUTS
    #    PROFILE-A list with as many elements as there are points to draw profiles for
    #       $PROFILE.DATA - A date x levels x variables matrix with atmospheric data for given point
    #       $LOCATION - A two element vector the lat/lon coordinates of the locations
    #       $FORECAST.DATE - Date and time of forecast

    if(!is.null(model.data$ensembles)) {
       if(length(unique(model.data$ensembles)) > 1) {
          stop("There appears to be more than one ensembles run in your data set!  BuildProfile can only handle one at a time.  Use SubsetNOMADS to split them up.")
       }
    }

    profile   <- NULL
 
    variables <- unique(model.data$variables)
    levels    <- unique(model.data$levels)
    times     <- unique(model.data$forecast.date)

    #Calculate distance from point of interest to each point on grid
    ll.all <- paste(model.data$lat, model.data$lon)
    ll.u  <- unique(ll.all)
    ll.tmp <- as.numeric(unlist(strsplit(ll.u, " ")))
    ll.arr <- t(array(ll.tmp, dim = c(2, length(ll.tmp)/2)))[, 2:1]

    #Standardize longitudes to permit model domain checks
    lon.tmp.point <- lon
    lon.tmp.model <- unique(model.data$lon)
    if(any(lon.tmp.model < 0)) {
        lon.tmp.point[lon.tmp.point < 0] <- lon.tmp.point[lon.tmp.point < 0] + 360
        lon.tmp.model[lon.tmp.model < 0] <- lon.tmp.model[lon.tmp.model < 0] + 360
    }
    lon.span <- range(lon.tmp.model)
    lat.tmp.point <- lat #Pedantic, but I'm human after all
    lat.span <- range(model.data$lat)
 
    profile.warning <- FALSE
    for(k in 1:length(lon)) { 
   
        if((lon.tmp.point[k] < min(lon.span)) | 
           (lon.tmp.point[k] > max(lon.span)) | 
           (lat.tmp.point[k] > max(lat.span)) | 
           (lat.tmp.point[k] < min(lat.span))) {
 
           warning(paste0("Profile point ", k, " of ", length(lon), " may be outside of the model domain, or it's simply close to (the antipode of) the prime meridian.\n", 
           "See GetProfile documentation for details."))
        }

        dist <- fields::rdist.earth(ll.arr, cbind(lon[k], lat[k]))
        profile.data <- array(NA,
            dim = c(length(levels), length(variables), length(times)))
    
        for(j in 1:length(times)) {
            if(spatial.average) {  
                s.i <- sort(dist, index = TRUE)$ix[1:points]
                sa.tmp <- as.numeric(unlist(strsplit(ll.u[s.i], " ")))
                sa.arr <- t(array(sa.tmp, dim = c(2, length(sa.tmp)/2)))[, 2:1]
                proj <- GEOmap::setPROJ(type = 2, LAT0 = lat[k], LON0 = lon[k])
                cart.pts <- GEOmap::GLOB.XY(sa.arr[, 2], sa.arr[, 1], proj)
        
                for(l in seq_len(length(levels))) {
                    for(m in seq_len(length(variables))) {
                        vl.i   <- which(model.data$levels == levels[l] & model.data$variables == variables[m] & model.data$forecast.date == times[j])
                        ll.vl  <- ll.all[vl.i]
                        val.vl <- model.data$value[vl.i]
                        d.i <- match(ll.u[s.i], ll.vl)             
                        layer.img <- cbind(cart.pts$x, cart.pts$y, val.vl[d.i])
                        if((all(layer.img[, 1] > 0) | all(layer.img[, 1] < 0)) | 
                           (all(layer.img[, 2] > 0) | all(layer.img[, 2] < 0)))
                        {
                            
                            profile.warning <- TRUE
                        }

                        profile.data[l, m, j] <- MBA::mba.points(layer.img, cbind(0, 0))[[1]][3]
                    }
                }
             } else { #Nearest grid node
                 d.i  <- which(ll.all == ll.u[which(dist == min(dist))] & model.data$forecast.date == times[j])
                 var.tmp <- model.data$variables[d.i]
                 lev.tmp <- model.data$levels[d.i]
                 val.tmp <- model.data$value[d.i]
                 for(l in 1:length(variables)) {
                     v.i <- which(var.tmp == variables[l])
                     profile.data[match(lev.tmp[v.i], levels), l, j]  <- as.numeric(val.tmp[v.i])
                 }
             }
         }

         profile[[k]] <- list(profile.data = profile.data, location = c(lon[k], lat[k]), forecast.date = times,
         variables = variables, levels = levels)
   } 

   if(profile.warning) {
       warning("A profile point fell outside of the model subset defined by the POINTS argument. 
           Increase POINTS to make sure you're interpolating rather than extrapolating.")
   }

   return(profile)
}

ModelGrid <- function(model.data, resolution, levels = NULL, variables = NULL, model.domain = NULL) {
    #Transform model data array into a grid with dimensions levels x variables x lon range x lat range
    #This should reduce the size of the returned data by removing redundant information
    #INPUTS
    #    MODEL.DATA - Data returned by ReadGrib
    #    RESOLUTION - Resolution of grid, in degrees if TYPE = "LATLON", in kilometers if TYPE = CARTESIAN, as a 2 element vector c(East-West, North-South)
    #    VARIABLES - variables to include in grid, if NULL, include all of them
    #    LEVELS - levels to include in grid, if NULL, include all of them
    #    MODEL.DOMAIN - vector c(LEFT LON, RIGHT LON, TOP LAT, BOTTOM LAT) of region to include in output. If NULL, include everything.
    #OUTPUTS
    #   FCST.GRID - A list with elements:
    #       $Z An array of dimensions levels x variables x lon x lat; each level x variable contains the model grid of data from that variable and level
    #       $X Vector of longitudes
    #       $Y Vector of latitudes
    #       $VARIABLES - the variables contained in the grid
    #       $LEVELS - the levels in the grid
    #       $MODEL.RUN.DATE - when the forecast model was run
    #       $FCST.DATE - what date the forecast is for
  
    if(length(resolution) != 2) {
        stop("\"resolution\" must be a 2 element vector: c(ZONAL RESOLUTION, MERIDIONAL RESOLUTION)")
    }
 
    #Check to make sure the grid is not too fine
    data.lat <- round(diff(model.data$lat), 11)
    data.lon <- round(diff(model.data$lon), 11)
    if(median(data.lon[data.lon != 0]) > resolution[1]) {
         warning(paste("The resolution you've chosen appears to be finer than the resolution of the model.  You chose", resolution[1], "but the model longitude list suggests that", median(data.lon[data.lon != 0]), "is the actual resolution.  This could cause unexpected behavior."))
    } 
    if(median(data.lat[data.lat != 0]) > resolution[2]) {
        warning(paste("The resolution you've chosen appears to be greater than the resolution of the model.  You chose", resolution[2], "but the model latitude list suggests that", median(data.lat[data.lat != 0]), "is the actual resolution.  This could cause unexpected behavior."))
    } 
    nodes.xy <- cbind(model.data$lon, model.data$lat)

    model.run.date <- unique(model.data$model.run.date)
    if(length(model.run.date) > 1) {
        stop("There appears to be more than one model run date in your data set! ModelGrid can only handle one at a time.")
    }

    if(!is.null(model.data$ensembles)) {
       if(length(unique(model.data$ensembles)) > 1) {
          stop("There appears to be more than one ensembles run in your data set!  ModelGrid can only handle one at a time.")
       }
    }

    fcst.date <- as.POSIXlt(unique(model.data$forecast.date), tz = "GMT")

    if(length(fcst.date) > 1) {
        stop("There appears to be more than one forecast in your data set! ModelGrid can only handle one at a time.")
    }

    if(is.null(variables)) {
        variables <- unique(model.data$variables) 
    }
  
    nomatch.ind <- is.na(match(variables, unique(model.data$variables)))
    if(sum(nomatch.ind) > 0) {
        warning(paste("The following variables are NOT present in the model data:", paste(variables[nomatch.ind], collapse = " ")))
        variables <- variables[!nomatch.ind]
    }

 
    if(is.null(levels)) {
        levels <- unique(model.data$levels)
    }

    nomatch.ind <- is.na(match(levels, unique(model.data$levels)))
    if(sum(nomatch.ind) > 0) {
        warning(paste("The following levels are NOT present in the model data:", paste(levels[nomatch.ind], collapse = " ")))
        levels <- levels[!nomatch.ind]
    }


    if(is.null(model.domain)) {
        model.domain <- c(min(nodes.xy[,1]), max(nodes.xy[,1]), max(nodes.xy[,2]), min(nodes.xy[,2]))
    }

    #Build grid

    grid <- list(x = seq(model.domain[1], model.domain[2], by = resolution[1]),
       y = seq(model.domain[4], model.domain[3], by = resolution[2]))
    
    fcst.grid <- list(z = array(rep(NA, length(grid$x) * length(grid$y) * length(variables) * length(levels)),
        dim = c(length(levels), length(variables), length(grid$x), length(grid$y))), 
        x = grid$x, y = grid$y, variables = variables, levels = levels, 
        model.run.date = model.run.date, fcst.date = fcst.date)

    #Put variables and levels into a series of layered images
    for(lvl in levels) {
        for(var in variables) {
             mi <- which(var == model.data$variables & lvl == model.data$levels) 
             if(length(mi) > 0) {
                 fcst.grid$z[which(lvl == fcst.grid$levels), which(var == fcst.grid$variables),,] <- fields::as.image(
                     as.numeric(model.data$value[mi]),
                     grid = grid,
                     x = nodes.xy[mi,])$z
              }
        }
    }

    return(fcst.grid)
}

MagnitudeAzimuth <- function(zonal.wind, meridional.wind) {
   #Given zonal (East-West) and meridional (North-South) wind speeds, calculate magnitude and azimuth.
   #INPUTS
   #    ZONAL.WIND - Wind East West, in meters per second, west negative
   #    MERIDIONAL.WIND - Wind North South, in meters per second, south negative
   #OUTPUTS
   #   MAGNITUDE - Wind magnitude, in meters per second
   #   AZIMUTH  - Wind azimuth, in degrees from north
   
   mag <- sqrt(zonal.wind^2 + meridional.wind^2)
   tmp.az <- (180/pi) * atan2(zonal.wind, meridional.wind)
   az <- tmp.az
   az[tmp.az < 0] <- 360 + tmp.az[tmp.az < 0]
   return(list(magnitude = mag, azimuth = az))
}

LinkExtractor <- function(url) {
    #Extracts links from a web page
    #INPUTS
    #    URL - Address of web page
    #OUTPUTS
    #    LINKS - A list of all the links on the page

    html.tmp <- readLines(url, warn = FALSE)

    hrefs    <- unlist(strsplit(html.tmp[which(grepl("href", html.tmp))],
       "</tr><tr><td>"))

    meat     <- stringr::str_extract(hrefs, "href=\".*\"")

    links <- stringr::str_replace_all(
        stringr::str_replace(meat, "href=", ""), 
        "\"", "")

    links <- stringr::str_replace_all(links, "http:", "https:")        

    links <- stringr::str_replace_all(links, ":80", ":443")

    return(links)
}

PlotWindProfile <- function(zonal.wind, meridional.wind, height, magnitude = NULL, magnitude.range = c(0, 50), 
    height.range = c(0, 50000), points = TRUE, lines = FALSE,
    radial.axis = TRUE, elev.circles = NULL, elev.labels = NULL, 
    radial.lines = NULL, colorbar = TRUE, colorbar.label = NULL, north.label = TRUE, invert = FALSE, ...) {
    #Plot a graphical representation of wind speed and direction, with lowest elevations closest to the center and highest elevations towards the edge.
    #INPUTS
    #    ZONAL.WIND - East-west wind speed.  Can be a vector or a list of vectors
    #    MERIDIONAL.WIND - North-South wind speed.  Can be a vector or a list of vectors.
    #    HEIGHT - Height above reference point.  Can be a vector or list of vectors.
    #    MAGNITUDE - If you want to plot different numeric values along the azimuths specified by ZONAL.WIND and MERIDIONAL.WIND, for example, if you want
    #       along-wind sound velocity.  Default NULL.
    #    MAGNITUDE.RANGE - Ranges of wind speed to plot, defaults to c(0, 50)
    #    HEIGHT.RANGE - Which heights to plot, defaults to c(0, 50000)
    #    POINTS - Do you want to plot data points as points?  You can use optional parameters to set the type and appearance of these points. Default TRUE.
    #    LINES - DO you want to plot your data as lines?  You can use optional parameters to set the type and appearance of these lines. Default FALSE.
    #    RADIAL.AXES - Plot an axis around the figure
    #    ELEV.CIRCLES - A vector of elevations to plot as a grid.  Default NULL
    #    ELEV.LABELS - Labels to put above each elevation circle
    #    RADIAL.LINES - A vector of azimuths to plot.  Default NULL
    #    COLORBAR - Whether to put a colorbar of wind magnitudes
    #    COLORBAR.LABEL - What to label your colorbar
    #    NORTH.LABEL - whether to label north. Default TRUE.
    #    INVERT - Whether to reverse the plot, so that higher elevations are towards the center and lower elevations are towards the outer edges.
    #OPTIONAL INPUTS (via ...)
    #    R.AXIS - radius of plot axis
    #    TICK.LEN - length of azimuth ticks
    #    R.AXIS.TICKS - Whether or not to put tick marks on the outer axis
    #    MAX.AZ - If plotting lines and the difference between two segments is greater than this value, interpolate between them to make things smooth
    #    COLOR.MAP - A list of colors to use, defaults to rainbow(n, 0, 5/6)
    #    N.COLS - Number of color bins in color map
    #    SUB.COL - Color of internal (elevation and azimuth) axes as a vector of length 2
    #    SUB.LTY - Type of internal axes, as a vector of length 2
    #    SUB.LWD - Width of internal axes, as a vector of length 2
    #    ELEV.LABELS.AZ - Which azimuth to plot elevation labels on. 
    #    POINT.CEX - size of points, if plotted
    #    PCH - Plot character of points, if plotted
    #    LTY - Line style, if lines are selected
    #    LWD - Line thickness, if lines are selected
    #    COLORBAR.TICK - Where to put labels on colorbar
 
    opts <- list(...)
    o.n <- names(opts)
    r <- 1

    if(! "r.axis" %in% o.n) {
        r.axis <- 1.2
    } else {
        r.axis <- opts$r.axis
    }

    if(! "tick.len" %in% o.n) {
        tick.len <- r / 20 
    } else {
        tick.len <- opts$tick.len
    }
    
    if(! "r.axis.ticks" %in% o.n) {
        r.axis.ticks <- TRUE
    } else {
        r.axis.ticks <- opts$r.axis.ticks
    }

    if(! "max.az" %in% o.n) {
        max.az <- 5
    } else {
        max.az <- opts$max.az
    }

    if(!"n.cols" %in% o.n) {
        n.cols <- 500
    } else {
        n.cols <- opts$n.cols
    }

    if(!"color.map" %in% o.n) {
       color.map <- rainbow(n.cols, start = 0, end = 5/6)           
    } else {
       color.map <- opts$color.map
    }

    if(!"sub.col" %in% o.n) {
        sub.col <- c("gray50", "gray60")
    } else {
       sub.col <- opts$sub.col
    }

    if(!"sub.lty" %in% o.n) {
        sub.lty <- c(2, 3)
    } else {
       sub.lty <- opts$sub.lty
    }

    if(!"sub.lwd" %in% o.n) {
        sub.lwd <- c(1, 1)
    } else {
        sub.lwd <- opts$sub.lwd
    }

    if(!"elev.labels.az" %in% o.n) {
        elev.labels.az <- 0
    } else {
        elev.labels.az <- opts$elev.labels.az
    }

    if(!"point.cex" %in% o.n) {
        point.cex <- 1
    } else {
        point.cex <- opts$point.cex
    }

    if(!"pch" %in% o.n) {
        pch <- 1
    } else {
       pch <- opts$pch
    }

    if(!"lty" %in% o.n) {
        lty <- 1
    } else {
        lty <- opts$lty
    }

    if(!"lwd" %in% o.n) {
        lwd <- 1
    } else {
       lwd <- opts$lwd
    }

    if(!"colorbar.tick" %in% o.n) {
        colorbar.tick <- seq(magnitude.range[1], magnitude.range[2], length.out = 6)
    } else {
        colorbar.tick <- opts$colorbar.tick
    }

    rng <- max(height.range) - min(height.range)

    plot(c(-r.axis - tick.len * 4, r.axis + tick.len * 4) + r/5, 
        c(-r.axis - tick.len * 4, r.axis + tick.len * 4), asp = 1, 
        type = "n", axes = FALSE, xlab = "", ylab = "") 

    if(radial.axis) {
        angs <- seq(0, 2 * pi, length.out = 1000)
        circ.axis <- list(x = r.axis * cos(angs), y = r.axis * sin(angs))
        lines(circ.axis)
        if(r.axis.ticks) {
            circ.ticks <- seq(0, 360, by = 45)[1:8] * pi/180
            tick.len <- r/20
            segments(r.axis * cos(circ.ticks), r.axis * sin(circ.ticks),
                r.axis * cos(circ.ticks) + tick.len * cos(circ.ticks),
                r.axis * sin(circ.ticks) + tick.len * sin(circ.ticks))
        }
   }

   if(north.label) { 
       text(0, r.axis + tick.len * 2, lab = "North")
   }

   interp.len <- 100  #Number of points in azimuth interpolation, if necessary

   if(! is.list(zonal.wind)) {
       zonal.wind <- list(zonal.wind)
       meridional.wind <- list(meridional.wind)
       height <- list(height)
   }

   color.scale <- seq(magnitude.range[1], magnitude.range[2], length.out = n.cols) #Wind color scale

   for(k in 1:length(zonal.wind)) {
       hgt.ind <- which(height[[k]] >= height.range[1] & height[[k]] <= height.range[2])
       maz <- MagnitudeAzimuth(zonal.wind[[k]][hgt.ind], meridional.wind[[k]][hgt.ind])
       hgt <- height[[k]][hgt.ind] - min(height.range)
       if(is.null(magnitude)) {
           mag <- maz$magnitude
       } else {
           mag <- magnitude[[k]][hgt.ind]
       }

       if(invert) {
           hgt <- rev(hgt)
       }
       if(points) {
           point.col <- c()
           if(is.null(magnitude)) {
               mag <- maz$magnitude
           } else {
               mag <- magnitude
           }
           for(j in 1:length(maz$azimuth)) {
                point.col <- append(point.col, color.map[which(abs(mag[j] - color.scale) == min(abs(mag[j] - color.scale)))])
           }
           points(r * hgt/rng * sin((pi / 180) * maz$azimuth),
               r * hgt/rng * cos((pi / 180) * maz$azimuth), 
               col = point.col, cex = point.cex, pch = pch)
       }
       if(lines) {
           m <- c()
           h <- c()
           azs <- c(maz$azimuth, tail(maz$azimuth, 1))
           tmp.az <- c()
           for(j in 1:length(maz$azimuth)) {
               if(abs(azs[j] - azs[j + 1]) < max.az) {
                   m <- append(m, mag[j])
                   h <- append(h, hgt[j])
                   tmp.az <- append(tmp.az, azs[j])
               } else { #Smooth out large azimuth swings (usually at low wind speeds)
                   #If we do not cross north
                   if((abs(azs[j] - (azs[j + 1] + 360)) > (abs(azs[j] - azs[j +1]))) & (abs(azs[j] + 360 - azs[j + 1]) > abs(azs[j] - azs[j + 1]))) {
                       tmp.az <- append(tmp.az, seq(azs[j], azs[j +1], length.out = interp.len))
                       #If we cross north, going clockwise
                   } else if (abs(azs[j] - (azs[j + 1] + 360)) < abs(azs[j] - azs[j + 1])) {
                       tmp.az <-  append(tmp.az, seq(azs[j], azs[j + 1] + 360, length.out = interp.len))
                   #If we cross north, going anticlockwise
                   } else {
                       tmp.az <-  append(tmp.az, seq(azs[j] + 360, azs[j + 1], length.out = interp.len))
                  }
                  m <- append(m, seq(mag[j], mag[j + 1], length.out = interp.len))
                  h <- append(h, seq(hgt[j], hgt[j + 1], length.out = interp.len))
               }
            }
            x <- r * h/rng * sin((pi / 180) * tmp.az)
            y <- r * h/rng * cos((pi / 180) * tmp.az)
            line.col <- c()
            for(j in 1:length(m)) {
               line.col <- append(line.col, color.map[which(abs(m[j] - color.scale) == min(abs(m[j] - color.scale)))])
            }
            segments(x[2:length(m)], y[2:length(m)], x[1:(length(m) - 1)], y[1:(length(m) - 1)], 
                col = line.col, lty = lty, lwd = lwd)
        }
   }

   #Altitude axes and labels
   text.loc <- elev.labels.az * (pi/180)
   if(!is.null(elev.circles)) {
       for(k in 1:length(elev.circles)) {
           h.axis <- list(x = r * ((elev.circles[k] - height.range[1])/rng) * cos(angs), y = r * ((elev.circles[k] - height.range[1])/rng) * sin(angs)) 
           lines(h.axis, lty = sub.lty[1], lwd = sub.lwd[1], col = sub.col[1])
           if(!is.null(elev.labels)) {
               text(r * sin(text.loc) * ((elev.circles[k] - height.range[1])/rng),  r * cos(text.loc) * ((elev.circles[k] - height.range[1])/rng), elev.labels[k], adj = c(0, 0))
           }
       }
   }

   #Radial axes and labels
   if(!is.null(radial.lines)) {
       segments(0, 0, r.axis * sin(radial.lines * (pi/180)), r.axis * cos(radial.lines * (pi/180)), lty = sub.lty[2], lwd = sub.lwd[2], col = sub.col[2])
   }
 
   #Make colorbar
   if(colorbar) {
       image(seq(r.axis + tick.len * 4, r.axis + tick.len * 5, length.out = 2),
           seq(-r.axis, r.axis, length.out = n.cols),
           array(seq(0, 10, length.out = n.cols),
           dim = c(1, n.cols)), add = TRUE, col = color.map) 
       scaling <- r.axis/(max(colorbar.tick) - min(colorbar.tick))
       text(rep(r.axis + tick.len * 4, length(colorbar.tick)), (colorbar.tick - min(colorbar.tick)) * 2 * scaling - r.axis,
           lab = colorbar.tick, pos = 4 )

       if(!is.null(colorbar.label)) {
           text(r.axis + tick.len * 3, 0, colorbar.label, srt = 90)
       }
    }
}     

SubsetNOMADS <- function(model.data, levels = NULL, variables = NULL, lon = NULL, lat = NULL,
ensembles = NULL, forecast.date = NULL, model.run.date = NULL) {
   
   #Subset the output of DODSGrab or ReadGrib to include only those parameters you are interested in
   #INPUTS
   #    MODEL.DATA     - Output of DODSGrab or ReadGrib
   #    LEVELS         - Which levels to keep, as vector
   #    VARIABLES      - Which variables to keep, as vector
   #    LON            - Which longitude points to keep, as vector, note this culd be very inefficient
   #    LAT            - Which latitude points to keep, as vector, note this culd be very inefficient
   #    ENSEMBLE       - Which ensembles runs to keep, as vector
   #    FORECAST.DATE  - Which forecast dates to keep, as vector.
   #    MODEL.RUN.DATE - Which model run dates to keep, as vector; not sure if this will ever be used 
   #OUTPUTS
   #    MODEL.DATA.SUB - Subsetted model.data per above criteria

   d.i <- rep(1, length(model.data$value))
 
   if(!is.null(levels)) { 
      d.i <- d.i & model.data$levels %in% levels 
   }

   if(!is.null(variables)) {
       d.i <- d.i & model.data$variables %in% variables
   }

   if(!is.null(lon)) {
       d.i <- d.i & model.data$lon %in% lon
   }

   if(!is.null(lat)) {
       d.i <- d.i & model.data$lat %in% lat
   }

   if(!is.null(ensembles)) {
       d.i <- d.i & model.data$ensembles %in% ensembles
   }

   if(!is.null(forecast.date)) {
       d.i <- d.i & as.numeric(model.data$forecast.date) %in% as.numeric(as.POSIXct(forecast.date))
   }

   if(!is.null(model.run.date)) {
       for(date.tmp in model.run.date) {
           d.i <- d.i & (date.tmp == model.run.date)
       }
   }

   list.names <- names(model.data)
   model.data.sub <- NULL

   for(list.name in list.names) {
      if(length(model.data[[list.name]]) == length(d.i)) {
          model.data.sub[[list.name]] <- model.data[[list.name]][d.i]
      } else {
          model.data.sub[[list.name]] <- model.data[[list.name]]
      }

   }
   return(model.data.sub)
}

Try the rNOMADS package in your browser

Any scripts or data that you put into this service are public.

rNOMADS documentation built on Sept. 6, 2022, 5:06 p.m.