R/forecast-readwrite.R

In coolTot/DtgRecup: Apirateur de contenus Web Hydro-météorologiques

##' Lecture et ecriture de grilles de previsions
##' 
##' \command{writeForecast} enregistre dans un fichier au format 
##' netCDF un jeu de prévisions pour une varaible donnée.
##' 
##' @param x un jeu de prévision de classe \sQuote{\var{Forecast}}. 
##' @param path le répertoire de suavegarde du jeu de préisions. 
##' 
<<<<<<< HEAD
##' @author Thomas Esclaffer <[email protected]@edf.fr>
##' @importFrom raster extent xres yres
##' @import ncdf4
=======
##' @import ncdf
##' @author Thomas Esclaffer <[email protected]@edf.fr>
>>>>>>> 53ae716ce8775c9fd36538db2ce00c6b918e4778
##' @export 
writeForecast <- 
  function(x, path=".") {
    ## Creation du nom de ficher
    file <- sprintf("%s_%s_%s.nc", toupper(x$model),
                    format(x$run,"%Y-%m-%d-%H", tz="UTC"), x$variable)
    file <- file.path(path, file)
    ## Préparation du fichier netCDF
<<<<<<< HEAD
    xt <- extent(x$data) 
    rx <- xres(x$data) 
    ry <- yres(x$data)
    vlon <- seq(from=xt[1]+rx/2, to=xt[2]-rx/2, by=rx)
    vlat <- seq(from=xt[3]+ry/2, to=xt[4]-ry/2, by=ry)    
    ## Description des dimensions
    lon <- ncdim_def("lon","degreesE", vlon)
    lat <- ncdim_def("lat","degreesN", vlat)
    ech <- ncdim_def("ech", "hours", x$ech)
    dim <- list(lat, lon, ech)
    ## Description des variables
    var <- ncvar_def(name = x$variable, units = x$units, 
                            dim = dim, missval=NA)
    ## Creation de l'archive netCDF
    ddim <- sapply(dim, with, len)
    if (any(dim(x$data) != ddim))
      stop("Dimensions des tableaux incompatibles !", dim(x$data),ddim)
    nc <- nc_create(file, var) 
    ncvar_put(nc, var, as.array(x$data))
    ## Attributs supplémentaires 
    ncatt_put(nc, 0, "model", x$model)
    ncatt_put(nc, 0, "run", 
                       format(x$run, "%Y%m%d%H", tz="UTC"))
    ##
    nc_close(nc)
=======
    xt <- raster::extent(x$data) 
    rx <- raster::xres(x$data) 
    ry <- raster::yres(x$data)
    vlon <- seq(from=xt[1]+rx/2, to=xt[2]-rx/2, by=rx)
    vlat <- seq(from=xt[3]+ry/2, to=xt[4]-ry/2, by=ry)    
    ## Description des dimensions
    lon <- ncdf::dim.def.ncdf("lon","degreesE", vlon)
    lat <- ncdf::dim.def.ncdf("lat","degreesE", vlat)
    ech <- ncdf::dim.def.ncdf("ech", "hours", x$ech)
    ## Description des variables
    var <- ncdf::var.def.ncdf(name = x$variable, units = x$units, 
                              dim = list(lon, lat, ech), missval=NA)
    ## Creation de l'archive netCDF
    ddim <- c(length(vlon), length(vlat), length(x$ech))
    if (any(dim(x$data) != ddim))
      stop("Dimensions des tableaux incompatibles !", dim(x$data),ddim)
    nc <- ncdf::create.ncdf(file, var) 
    ncdf::put.var.ncdf(nc, var, as.array(x$data))
    ## Attributs supplémentaires 
    ncdf::att.put.ncdf(nc, 0, "model", x$model)
    ncdf::att.put.ncdf(nc, 0, "run", 
                       format(x$run, "%Y%m%d%H", tz="UTC"))
    ##
    ncdf::close.ncdf(nc)
>>>>>>> 53ae716ce8775c9fd36538db2ce00c6b918e4778
    invisible(file)
  }

##' @rdname writeForecast
##' @description 
##' \command{readForecast} charge le jeu de prévision contenu 
##' dans un fichier netCDF, produit à l'aide de 
##' \command{writeForecast}.
##' 
##' @param file le nom du fichir de prévisions à charger. 
##' 
<<<<<<< HEAD
##' @import ncdf4
##' @importFrom raster brick
##' 
=======
>>>>>>> 53ae716ce8775c9fd36538db2ce00c6b918e4778
##' @export
readForecast <- 
  function(file) 
  {
<<<<<<< HEAD
    if (grepl("^ftp://", file)) {
      tmp <- file.path(tempdir(), basename(file))
      download.file(file, tmp, mode = "wb", quiet = TRUE)
      file <- tmp
    }
    if (!file.exists(file)) stop("Fichier inexistant !")
    nc <- nc_open(file, write=FALSE)
=======
    if (!file.exists(file)) stop("Fichier inexistant !")
    nc <- ncdf::open.ncdf(file, write=FALSE)
>>>>>>> 53ae716ce8775c9fd36538db2ce00c6b918e4778
    ## Recupération des propriétés des variables
    var <- sapply(nc$var, with, name)
    units <- sapply(nc$var, with, units)
    ech <- nc$dim$ech$vals
    ## Récupération des attributs
<<<<<<< HEAD
    model <- ncatt_get(nc,0,"model")$value
    run <- ncatt_get(nc,0,"run")$value
    run <- as.POSIXct(run,"%Y%M%d%H", tz="UTC")
    ## Systeme decoordonnées 
    if (nc$dim[[1]]$name == "lat") {
      crs <- CRS("+init=epsg:4326")
    } else crs <- CRS("+init=epsg:27572")
    ## Mise en forme des variables
    data <- ncvar_get(nc)
    yres<- diff(nc$dim[[1]]$vals)[1]
    xres<- diff(nc$dim[[2]]$vals)[1]
    rd <- brick(data, 
                xmn=min(nc$dim[[2]]$vals)-xres/2, 
                xmx=max(nc$dim[[2]]$vals)+xres/2,
                ymn=min(nc$dim[[1]]$vals)-yres/2, 
                ymx=max(nc$dim[[1]]$vals)+yres/2,
                crs=crs)
=======
    model <- ncdf::att.get.ncdf(nc,0,"model")$value
    run <- ncdf::att.get.ncdf(nc,0,"run")$value
    run <- as.POSIXct(run,"%Y%M%d%H", tz="UTC")
    ## Systeme decoordonnées 
    if (nc$dim[[1]]$name == "lon") {
      crs <- sp::CRS("+init=epsg:4326")
    } else crs <- sp::CRS("+init=epsg:27572")
    ## Mise en forme des variables
    data <- ncdf::get.var.ncdf(nc)
    xres<- diff(nc$dim[[1]]$vals)[1]
    yres<- diff(nc$dim[[2]]$vals)[1]
    rd <- raster::brick(data, 
                        xmn=min(nc$dim[[1]]$vals)-xres/2, 
                        xmx=max(nc$dim[[1]]$vals)+xres/2,
                        ymn=min(nc$dim[[2]]$vals)-yres/2, 
                        ymx=max(nc$dim[[2]]$vals)+yres/2,
                        crs=crs)
>>>>>>> 53ae716ce8775c9fd36538db2ce00c6b918e4778
    ## 
    forecast(data=rd, model = model, run = run, variable = var,
             echeances = ech, units=units)
  }