Create raster from SAFRAN csv file

safran2raster.R

# safran 2 raster

library(data.table)
library(raster)
library(sf)


# ##### FONCTIONS #####

# genere 1 raster à partir de DT(data.table) pour la variable et la date spécifiées
DT2rast <- function(extr_date, DT, varname) {
  varlist <- c("LAMBX", "LAMBY", varname)
  DTi <- DT[DATE==extr_date, ..varlist]
  dfi <- as.data.frame(DTi)
  dfi$LAMBX <- dfi$LAMBX * 100
  dfi$LAMBY <- dfi$LAMBY * 100
  rasterFromXYZ(as.data.frame(dfi), crs=CRS("+init=EPSG:27572"), digits=0)
}

# genere 1 stack raster à partir de DT(data.table)
# pour la variable, la plage temporelle et l'étendue spécifiées
safran2rast <- function(DT, varname, d_date, f_date, etendue_l2e=NULL) {
  varlist <- c("LAMBX", "LAMBY", "DATE", varname)
  seq_dates <- strftime(seq(d_date, f_date, by="days"), "%Y%m%d")
  if (!is.null(etendue_l2e)) {
    # subset etendue Lambert2 en hm
    ethm <- round(etendue_l2e/100, 0)
    DT <- DT[LAMBX>ethm[1] & LAMBX<ethm[3] & LAMBY>ethm[2] & LAMBY<ethm[4]]
    DT <- DT[DATE %in% seq_dates,..varlist]
  } else {
    DT <- DT[DATE %in% seq_dates,..varlist]
  }
  list_raster <- lapply(seq_dates, DT2rast, DT, varname)
  stack(list_raster)
}

# genere grille points à partie de DT(data.table) pour la date spécifiée
DT2sp <- function(extr_date, DT) {
  DTi <- DT[DATE==extr_date]
  pts <- as.data.frame(DTi)
  pts$coordx <- pts$LAMBX * 100
  pts$coordy <- pts$LAMBY * 100
  geomdef <- paste0("POINT(", pts$coordx, " ", pts$coordy,")") 
  sfc <- st_as_sfc(geomdef, crs=27572)
  st_geometry(pts) <- sfc
  pts
}


# ##### EXEMPLES D'UTILISATION #####

# ---- parametre : fichier en entree ----

setwd("D:/MeteoFrance/Safran-Isba données quotidiennes")
f_safran <- "SIM2_2010_2017.csv"

DT_safran <- fread(f_safran, header = TRUE)

# etendue de l'Occitanie en Lambert 2 etendu ?
sf_reg <- st_read("regions-20180101-shp/regions-20180101.shp")
sf_occ <- sf_reg[sf_reg$nom=="Occitanie",]
sf_occ_l2 <- st_transform(sf_occ, 27572)
# etendue format xmin, ymin, xmax, ymax
etendue <- st_bbox(sf_occ_l2)
# dates debut et fin
jour1 <- as.Date("2014-01-01")
jourN <- as.Date("2014-12-31")


# ---- exemple 1 : process stack (multitemporel, 1 variable, etendue region) ----

stack_tm <- safran2rast(DT_safran, varname="T_Q", d_date=jour1, f_date=jourN, 
                         etendue_l2e=etendue)
stack_tn <- safran2rast(DT_safran, varname="TINF_H_Q", d_date=jour1, f_date=jourN, 
                         etendue_l2e=etendue)
stack_tx <- safran2rast(DT_safran, varname="TSUP_H_Q", d_date=jour1, f_date=jourN, 
                         etendue_l2e=etendue)
stack_rr <- safran2rast(DT_safran, varname="PRELIQ_Q", d_date=jour1, f_date=jourN, 
                         etendue_l2e=etendue)

writeRaster(stack_tm, "output_safran/tm_occitanie_2014.tif", overwrite=T)
writeRaster(stack_tn, "output_safran/tn_occitanie_2014.tif", overwrite=T)
writeRaster(stack_tx, "output_safran/tx_occitanie_2014.tif", overwrite=T)
writeRaster(stack_rr, "output_safran/rr_occitanie_2014.tif", overwrite=T)

# ... extraire temp + prec montpellier
cefe_l2e <- matrix(c(723548, 1849625), ncol=2)
cefe_stats <- data.frame(
  dates_ = strftime(seq(jour1, jourN, by="days"), "%Y%m%d"),
  tm = as.numeric(extract(stack_tm, cefe_l2e)),
  tn = as.numeric(extract(stack_tn, cefe_l2e)),
  tx = as.numeric(extract(stack_tx, cefe_l2e)),
  rr = as.numeric(extract(stack_rr, cefe_l2e))
)
write.csv2(cefe_stats, "output/cefe_stats.csv", row.names = F)


# ---- exemple 2 : process raster (1 jour, 1 variable, tte la France) ----

raster_tm_0703 <- DT2rast("20140929", DT_safran, "T_Q")
raster_tn_0703 <- DT2rast("20140929", DT_safran, "TINF_H_Q")
raster_tx_0703 <- DT2rast("20140929", DT_safran, "TSUP_H_Q")
raster_rr_0703 <- DT2rast("20140929", DT_safran, "PRELIQ_Q")

writeRaster(raster_tm_0703, "output_safran/tm_france_20140929.tif")
writeRaster(raster_tn_0703, "output_safran/tn_france_20140929.tif")
writeRaster(raster_tx_0703, "output_safran/tx_france_20140929.tif")
writeRaster(raster_rr_0703, "output_safran/rr_france_20140929.tif")


# ---- exemple 3 : process grille points (1 jour, N variables, tte la France) ----

pts_20180201 <- DT2sp("20140929", DT_safran)
st_write(pts_20180201, "output_safran/data_20140929.shp", delete_layer = T)