R/ncdf_jra55_minmax_surf.R

#' Download and extract the Japanese Re-Analysis 55 Years (JRA-55) minimum and maximum temperature at 2 meters
#' 
#' This function download JRA-55 minimum and maximum temperature at 2 meters above ground and extract it 
#' over a specified region.
#' 
#' @param start_date Start date to be downloaded, a character in the form "YYYY-MM-DD".
#' @param end_date End date to be downloaded, a character in the form "YYYY-MM-DD".
#' @param output.dir Full path to the directory where the downloaded and extracted files will be saved;
#' default to the working directory.
#' @param grads.bin Full path name of GrADS executable. For Mac and Linux, you can use \code{"grads"}
#' if GrADS is in your PATH environment variable.
#' @param range.lon A numeric vector of length 2; minimum and maximum longitude of the area to be extracted;
#' default c(-180, 180).
#' @param range.lat A numeric vector of length 2; minimum and maximum latitude of the area to be extracted;
#' default c(-89.570, 89.570).
#' @param tstep The time steps of the extracted data. Useful values are \code{"daily"} (aggregated to daily data)
#' and \code{"hourly3"} (3-hourly data).
#' @param username Your user ID for JMA Data Dissemination System (JDDS).
#' @param password Your password for JDDS.
#' 
#' @return A directory named "GRIB" and "Extracted" are created. "GRIB" contains the downloaded GRIB data and
#' "Extracted" contains the extracted data in NetCDF format over the area that you specified.
#' 
#' @examples
#' 
#' \dontrun{
#' temp_jra55_minmax_surf(
#'         start_date = "2019-05-26",
#'         end_date = "2019-05-26",
#'         output.dir = "~/Desktop/DATA/JRA55",
#'         grads.bin = "grads",
#'         range.lon = c(42, 52),
#'         range.lat = c(-26, -11.5),
#'         tstep = "daily",
#'         username = "my.UID",
#'         password = "my.PWD"
#'      )
#' }
#' 
#' @export

temp_jra55_minmax_surf <- function(
                                start_date = "2019-05-01",
                                end_date = "2019-05-31",
                                output.dir = getwd(),
                                grads.bin = "grads",
                                range.lon = c(-180, 180),
                                range.lat = c(-89.570, 89.570),
                                tstep = c("daily", "hourly3"),
                                username = NA, password = NA
                            )
{
    download_jra55_minmax_surf(start_date, end_date,
                               output.dir,
                               username, password)

    extarct_jra55_minmax_surf(start_date, end_date,
                              output.dir, grads.bin,
                              range.lon, range.lat,
                              tstep)
}

#' Extract the Japanese Re-Analysis 55 Years (JRA-55) minimum and maximum temperature at 2 meters
#' 
#' This function extract an existing downloaded GRIB data from JRA-55 minimum and maximum temperature at 2 meters above ground 
#' over a specified region.
#' 
#' @param start_date Start date to be extracted, a character in the form "YYYY-MM-DD".
#' @param end_date End date to be extracted, a character in the form "YYYY-MM-DD".
#' @param output.dir Full path to the directory where the extracted files will be saved.
#' @param grads.bin Full path name of GrADS executable. For Mac and Linux, you can use \code{"grads"}
#' if GrADS is in your PATH environment variable.
#' @param range.lon A numeric vector of length 2; minimum and maximum longitude of the area to be extracted.
#' @param range.lat A numeric vector of length 2; minimum and maximum latitude of the area to be extracted.
#' @param tstep The time steps of the extracted data. Useful values are \code{"daily"} (aggregated to daily data)
#' and \code{"hourly3"} (3-hourly data).
#' 
#' @return The extracted data are located under the directory called "Extracted".
#' 
#' @examples
#' 
#' \dontrun{
#' extarct_jra55_minmax_surf(
#'         start_date = "2019-05-26",
#'         end_date = "2019-05-26",
#'         output.dir = "~/Desktop/DATA/JRA55",
#'         grads.bin = "grads",
#'         range.lon = c(42, 52),
#'         range.lat = c(-26, -11.5),
#'         tstep = "daily"
#'      )
#' }
#' 
#' @export

extarct_jra55_minmax_surf <- function(
                                    start_date, end_date,
                                    output.dir, grads.bin,
                                    range.lon, range.lat, tstep
                                )
{
    if(WindowsOS()){
        if(!file.exists(grads.bin))
            stop(paste(grads.bin, "not found!"))
    }else{
        pthgs <- split_path(grads.bin)
        if(length(pthgs) == 1){
            ret <- system(paste('which', grads.bin), intern = TRUE)
            if(length(ret) == 0) stop("grads not found!")
        }else{
            if(!file.exists(grads.bin))
                stop(paste(grads.bin, "not found!"))
        }
    }

    grads_wraper_jra55_minmax_surf(start_date, end_date, output.dir, grads.bin)
    
    tstep <- tstep[1]
    hours <- time_step_3hour(start_date, end_date)
    dir.month <- format(hours, "%Y%m")

    ret <- lapply(unique(dir.month), function(mo){
        dir.mon <- file.path(output.dir, "GRIB", mo)
        get_jra55_minmax_surf(dir.mon, range.lon, range.lat, output.dir, tstep)
    })

    invisible()
}

get_jra55_minmax_surf <- function(
                            dir.month, range.lon, range.lat,
                             dir.out, tstep
                        )
{
    dir.ext <- file.path(dir.out, "Extracted", tstep)
    if(!dir.exists(dir.ext)) dir.create(dir.ext, recursive = TRUE)

    file.tmax <- file.path(dir.month, "tmax.nc")
    if(!file.exists(file.tmax))
        stop(paste(file.tmax, "not found!"))
    create_nc_jra55_minmax_surf("tmax", file.tmax, dir.ext, range.lon, range.lat, tstep)

    file.tmin <- file.path(dir.month, "tmin.nc")
    if(!file.exists(file.tmin))
        stop(paste(file.tmin, "not found!"))
    create_nc_jra55_minmax_surf("tmin", file.tmin, dir.ext, range.lon, range.lat, tstep)

    invisible()
}

create_nc_jra55_minmax_surf <- function(
                                    varid, ncfile, ncdir,
                                    range.lon, range.lat, tstep
                                )
{
    temp <- read_jra55_minmax_surf(ncfile, range.lon, range.lat, tstep)
    dx <- ncdim_def("Longitude", "degrees_east", temp$x)
    dy <- ncdim_def("Latitude", "degrees_north", temp$y)
    vars <- switch(varid,
            "tmax" = list("maximum temperature", max),
            "tmin" = list("maximum temperature", min)
        )
    longname <- switch(tstep,
            "daily" = paste("daily", vars[[1]], "computed from 3 hourly JRA55"),
            "hourly3" = paste("3-hourly", vars[[1]], "JRA55")
        )
    outnc.z <- ncvar_def(varid, "degC", list(dx, dy), -9999, longname, "float", compression = 9)

    dir.temp <- file.path(ncdir, varid)
    if(!dir.exists(dir.temp)) dir.create(dir.temp, recursive = TRUE)

    if(tstep == "daily"){
        idx <- split(seq_along(temp$t), format(temp$t, "%Y%m%d"))
        jour <- names(idx)
        ret <- lapply(seq_along(idx), function(j){
            fileout <- file.path(dir.temp, paste0(varid, "_", jour[j], ".nc"))
            tmp <- temp$z[, , idx[[j]]]
            tmp <- apply(tmp, 1:2, vars[[2]], na.rm = TRUE)
            tmp[is.na(tmp) | is.nan(tmp)] <- -9999
            nc <- nc_create(fileout, outnc.z)
            ncvar_put(nc, outnc.z, tmp)
            nc_close(nc)
        })
    }

    if(tstep == "hourly3"){
        tm <- format(temp$t, "%Y%m%d%H")
        ret <- lapply(seq_along(tm), function(j){
            fileout <- file.path(dir.temp, paste0(varid, "_", tm[j], ".nc"))
            tmp <- temp$z[, , j]
            tmp[is.na(tmp) | is.nan(tmp)] <- -9999
            nc <- nc_create(fileout, outnc.z)
            ncvar_put(nc, outnc.z, tmp)
            nc_close(nc)
        })
    }

    invisible()
}

read_jra55_minmax_surf <- function(ncfile, range.lon, range.lat, tstep)
{
    nc <- nc_open(ncfile)
    x <- nc$dim[[1]]$vals
    y <- nc$dim[[2]]$vals
    tm <- nc$dim[[4]]$vals
    units <- nc$dim[[4]]$units
    z <- ncvar_get(nc, varid = "var")
    nc_close(nc)
    x <- ((x + 180) %% 360) - 180
    orig <- paste0(strsplit(units, " ")[[1]][3:4], collapse = " ")
    orig <- as.POSIXct(orig, tz = "GMT", format = "%Y-%m-%d %H:%M")
    tm <- orig + tm * 60

    ix <- x >= range.lon[1] & x <= range.lon[2]
    iy <- y >= range.lat[1] & y <= range.lat[2]
    x <- x[ix]
    y <- y[iy]
    z <- z[ix, iy, ] - 273.15
    ox <- order(x)
    x <- x[ox]
    z <- z[ox, , ]
    i2 <- duplicated(x)
    if(any(i2)){
        x <- x[!i2]
        z <- z[!i2, , ]
    }

    list(t = tm, x = x, y = y, z = z)
}
rijaf-iri/CDTDownload documentation built on June 5, 2019, 12:37 a.m.