Here is my code
setwd("C:/CARTE_BAD")
# Charger les libraries --------------------------------------------------
#require(pacman)
#pacman::p_load(terra, sf, prettymapr,raster,geodata,RSAGA,tidyverse,
#gstat,glue,fs,readxl,RColorBrewer,tmap,rgdal,sp,pals,
#viridis)
library(terra) # Analyse données spatiales avec des données vectorielles
#(points, lignes, polygones) et raster (grille).
library(sf) # dédié à la manipulation, la transformation et
#l'analyse de données spatiales
library(prettymapr) # package pour les echelles et orientations
library(raster) # Manipulation des données raster
library(geodata) # pour le téléchargement de données géographiques
# à utiliser dans l'analyse spatiale et la cartographie
library(RSAGA) # analyse de terrain du système d'information géographique (SIG)
#"SAGA" (Système d'analyses géoscientifiques automatisées)
library(tidyverse) # Manipulation des données
library(gstat) # Packages pour les statistiques spatiales, modelisation
#geostatistique
library(fs) # Manipulation de donnees
library(readxl) # lire les fichiers excel
library(RColorBrewer)# Package pour les couleurs
library(pals) # Package pour les couleurs
library(viridis) # package pour les couleurs
library(tmap) # visualisation des données
library (tmaptools)
library(leaflet)
library(dplyr)
library(mapview)
# Creer le dossier de sortie des cartes -----------------------------------
rm(list = ls())
options(scipen = 999)
#dir.create("carte_du_bulletin/Aout-3-2022")
# Charger le fichier excel -----------------------------------------------
#tble <- read_excel('Postes1.xls')
#tble
#colnames(tble)
tble <- read_excel('RESULTAT_2022.xls',sheet = 24)
tble
colnames(tble)
# Change the colnames
cols <- c('PLUIE_DECA', 'CUML_PLUIE_DECA','DIFF_PLUI_DEC_N',
'DIFF_CUML_PLUIE_N','RU60mm','BHC_DEC')
cols
tble <- tble[,c(1,2,3, grep(paste0(cols, collapse = '|'), colnames(tble)))]
tble
# Telecharger le contour CIV ---------------------------------------------
civs <- geodata::gadm(country = 'CIV', level = 0, path = 'tmpr')
plot(civs, border = 'blue')
points(tble$LONG, tble$LAT, pch = 16, col = 'red')
# To project the shapefile ------------------------------------------------
civs <- terra::project(civs, '+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs')
ci_shp=st_read("ci_contour.shp")
ci_shp <-st_transform(ci_shp ,crs = st_crs("+proj=longlat +datum=WGS84 +no_defs"))
plot(ci_shp)
ocean_shp <- st_read("ocean.shp")
ocean_shp <-st_transform(ocean_shp ,crs = st_crs("+proj=longlat +datum=WGS84 +no_defs"))
plot(ocean_shp)
# Convert the table to shapefile -------------------------------------------
pnts <- tble
pnts <- drop_na(pnts)
#Remplacer les valeurs manquantes par -99
#pnts[is.na(pnts)] <- -9999
pnts1 <- st_as_sf(pnts, coords = c('LONG', 'LAT'), crs = st_crs(4326))
pnts <- st_transform(pnts1, st_crs(32630))
pnts <- terra::vect(pnts)
#tmap_mde('view')
tm_add
map1 <-tm_shape(ci_shp)+
tm_polygons(alpha = 0)+
tm_graticules()+
tm_shape(pnts1)+
#tm_bubbles(size = "PLUIE_DECA", col = "red", border.col = "black",
#title.size = "Precipitation \n décadaire [mm]",
#alpha = 0.7)
tm_dots(col="PLUIE_DECA", palette = "RdBu", stretch.palette = FALSE,
title="Precipitation \n décadaire [mm]", size=0.2) +
tm_text("PLUIE_DECA", col = "darkgreen", just="top", xmod=.5, ymod = 0.7, size = 0.9)+
tm_legend(legend.outside=TRUE)+
#tm_compass(position = c("right", "top"))+
tm_scale_bar()
map1
# Lecture et importation des stations -------------------------------------
Stations <- read_excel("Stations.xlsx")
crs_contour <- CRS(SRS_string="OGC:CRS84")
Stations$Long <- as.numeric(Stations$Long)
Stations$Lat <- as.numeric(Stations$Lat)
coordinates(Stations) <- ~Long + Lat
proj4string(Stations) <- crs_contour
class(Stations)
head(Stations, 10)
STATIONS_CRS<-spTransform(Stations,CRS(SRS_string="OGC:CRS84"))
STATIONS_CRS <-sf::st_as_sf(STATIONS_CRS)
# Make IDW ----------------------------------------------------------------
x.range <- terra::ext(civs)[1:2]
y.range <- terra::ext(civs)[3:4]
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 5000),
y = seq(from = y.range[1], to = y.range[2], by = 5000))
coordinates(grd) <- ~ x + y
gridded(grd) <- TRUE
raster::crs(grd) <- '+proj=utm +zone=30 +datum=WGS84 +units=m +no_defs'
pnts <- as(pnts, 'Spatial')
colnames(pnts@data) <- c('stt', 'pluie_deca', 'pluie_cumul','Dff_pluie',
'Dff_cumul','RU60mm','BHC')
head(pnts)
#-----------------------------------------
# PLUIE DECADAIRE
#-----------------------------------------
idw.d_pd <- gstat::idw(pluie_deca ~ 1, pnts, grd)
idw.d_pd <- raster::raster(idw.d_pd)
idw.d_pd <- rast(idw.d_pd)
idw.d_pd <- terra::crop(idw.d_pd, civs) %>% terra::mask(., civs)
idw.d_pd <- terra::project(idw.d_pd, '+proj=longlat +datum=WGS84 +no_defs')
#dir.create('raster')
terra::writeRaster(idw.d_pd, 'raster/idw_pd.tif')
# To download -------------------------------------------------------------
# SRTM --------------------------------------------------------------------
srtm <- geodata::elevation_30s(country = 'CIV', path = 'tmpr')
srtm <- terra::project(srtm, crs(idw.d_pd))
plot(srtm)
terra::writeRaster(srtm, 'raster/srtm.tif', overwrite = T)
# GWR ---------------------------------------------------------------------
env <- rsaga.env(path = 'C:/SAGA/saga-9.0.1_x64')
fle.srt <- 'raster/srtm.tif'
fle.inp <- 'raster/idw_pd.tif'
fle.out <- 'raster/gwr_pd.tif'
rsl <- rsaga.geoprocessor(
lib = 'statistics_regression',
module = 'GWR for Grid Downscaling',
param = list(PREDICTORS = fle.srt,
REGRESSION = fle.out,
DEPENDENT = fle.inp),
env = env
)
rst <- terra::rast(fle.out)
plot(rst,col = rainbow(25))
#m <- c(-1,0.25, 0.3, 0.4, 0.5, 1)
#vegc <- classify(ndvi, m)
#vegc_c <- crop(vegc, Cuenca_t)
# CODE COULEUR ------------------------------------------------------------
# Définir les couleurs pour chaque intervalle de valeurs
#colors <- c("#a50026","#d73027","#f46d43","#fdae61","#fee090","#ffffbf",
#"#e0f3f8","#abd9e9","#74add1","#4575b4","#313695")
# Définir les valeurs de pluviométrie et les couleurs correspondantes
VALEURS <- c(0, 10, 20, 30, 40, 50, 60, 75, 100, 150, 200, 250, 300)
# CARTE ---------------------------------------------------
#rev(topo.colors(256))
pluie_decadaire <- tm_shape(rst) +
tm_raster(col = "gwr_pd.tif", style = 'fixed',
n = 12, title = "Pluie décadaire\n[en mm]",
palette = "Spectral",
breaks = VALEURS,
#labels = c(seq(0, 600, 100), "> 700"),
midpoint = FALSE,
legend.reverse = TRUE, alpha = 1.0) +
tm_shape(STATIONS_CRS) +
tm_dots(col = "black", size = 0.1) +
tm_text("Stations", col = "black", size = 0.6, just = "top", ymod = 0.75) +
tm_shape(ci_shp) +
tm_borders(col = "black", lwd = 0.7) +
tm_shape(ocean_shp)+
tm_fill("lightskyblue1")+
tm_borders(col = "black", lwd = 0.7) +
tm_graticules(alpha = 0.4, lwd = 0.5, labels.size = 0.5) +
tm_compass(type = "8star",
position = c("RIGHT", "TOP"),
show.labels = 2,
text.size = 0.35) +
tm_scale_bar(position = c(0.4,0)) +
tm_logo("C:/CARTE_BAD/Image1.png",
height = 1.5,
halign = "bottom",
margin = 0.5,
position = c(0,0.9),
just = NA) +
tm_logo("C:/CARTE_BAD/ISO9001.PNG",
height = 1.8,
halign = "bottom",
margin = 0.1,
position = c(0,.1),
just = NA) +
tm_layout(
legend.format = list(text.separator = "-"),
main.title.fontface = 2,
main.title.position = 0.1,
main.title.size = 0.9,
panel.labels = c("Précipitation décadaire"),
panel.label.color = "darkslateblue",
inner.margins=c(0,0,0,0)) +
tm_xlab("Longitude (°W)") +
tm_ylab("Latitude (°N)") +
# Specify the variable for faceting
tm_facets(free.scales = TRUE) +
# Custom legend with fixed labels and colors
tm_legend(
outside = TRUE,
hist.width = 2,
legend.position = c(0.10, 0.2),
legend.bg.alpha = 1,
title.size = .9,
text.size = 0.7,
legend.bg.color = "gray90",
legend.frame = "gray90")
pluie_decadaire
tmap_save(pluie_decadaire, "C:/CARTE_BAD/carte_du_bulletin/pluie_deca_3_aout_22.png",
width = 25, height = 15, units = "cm")