R/rasterGFS.R

In meteoForecast: Numerical Weather Predictions

rasterGFS <- function(var,
                     day=Sys.Date(), run='00',
                     frames='complete',
                     box = NULL,
                     resolution = NULL,
                     names = NULL,
                     remote=TRUE, use00H = FALSE, ...) {

    ## Model initialization time
    run <- match.arg(run, mfRuns('gfs'))
    ## Time Frames
    if (remote) {
        ## MeteoGalicia implements netCDF Time Subset. Therefore,
        ## `frames` is a character.
        if (frames == 'complete') frames <- '&temporal=all'
        else {
            present <- as.POSIXct(paste0(as.character(day),
                                         as.numeric(run),
                                         ':00:00Z'))
            ff <- present
            lf <- present + (as.integer(frames) - 1) * 3600
            frames <- paste0('&time_start=',
                             format(ff,
                                    '%Y-%m-%dT%H:%M:%SZ'),
                             '&time_end=',
                             format(lf,
                                    '%Y-%m-%dT%H:%M:%SZ')
                             )
        }
    } else { ## remote=FALSE
        maxFrames <- mfHorizon("gfs") / mfTRes("gfs") + 1
        if (frames == 'complete') {
            frames <- seq(1, maxFrames, 1)
        } else {
            frames <- seq(1, min(frames, maxFrames), by=1)
        }
    }
        
    ## Name of files to be read/stored
    ncFile <- paste0(paste(var, ymd(day), run, 
                           sep='_'), '.nc')
        
    if (remote) {
        ## Meteogalicia provides a multilayer
        ## file with all the time frames
        completeURL <- composeURL(var, day, run,
                                  box, frames, resolution,
                                  service = 'gfs')
        message('Downloading data from ', completeURL)
        success <- try(download.file(completeURL, ncFile,
                                     mode='wb', quiet = TRUE),
                       silent=TRUE)
        if (inherits(success, "try-error")) {
            stop('Data not found. Check the date and variables name')
        } else { ## Download Successful!
            message('File(s) available at ', tempdir())
        } ## End of Remote
    } else {}
    ## Read files
    suppressWarnings(capture.output(bNC <- stack(ncFile)))
    ## Use frames with local files
    if (remote==FALSE) bNC <- bNC[[frames]]
    ## Convert into a RasterBrick
    b <- brick(bNC)
    ## Get values in memory to avoid problems with time index and
    ## projection
    b[] <- getValues(bNC)
    ## Projection parameters are either not well defined in the
    ## NetCDF files or incorrectly read by raster.
    projection(b) <- mfProj4('gfs')

    ## Use box specification with local files
    if (!is.null(box) & remote==FALSE){
        if (requireNamespace('sf', quietly=TRUE)) {
            extPol <- as(extent(box), 'SpatialPolygons')
            proj4string(extPol) <- '+proj=longlat +ellps=WGS84'
            extPol <- as(sf::st_transform(sf::st_as_sf(extPol),
                                          sf::st_crs(projection(b))),
                         "Spatial")
            b <- crop(b, extent(extPol))
        } else {
            warning("you need package 'sf' to use 'box' with local files")
        }
    }
    ## Time index
    hours <- seq(0, nlayers(b) - 1) * 3600 * mfTRes("gfs")
    tt <- hours + as.numeric(run)*3600 + as.POSIXct(day, tz='UTC')
    attr(tt, 'tzone') <- 'UTC'
    b <- setZ(b, tt)
    ## Names
    if (is.null(names)) names(b) <- format(tt, 'd%Y-%m-%d.h%H')
    ## Here it goes!
    b
}