#
# Géotraitments vectoriels avec R
# 
# Auteur  : Philippe Lejeune
# Version : 2022/09/13
#
#=============================================================


rm(list=ls()) # réinitialisation de la mémoire

# 1. Chargement des librairies --------------------

library(sf)
library(terra)
library(dplyr)
library(tidyverse) # mise en forme de tableau
library(lwgeom) # correction de géométrie
library(nngeo) # tracé d'ellipses, matrice de distance
library(DescTools) # applique 1 rotation à des objets


# 2. Lecture/écriture et principales propriétés des couches vectorielles ----------------



# !!!! utiliser "/" à la place de  "\" !!!!
path0="C:/PL/01_COURS/tthr_2022/R01_2022"
path_in=paste0(path0,"/data_in")
path_out=paste0(path0,"/output")

# Lire une couche vectorielle (-> objet sf)
f_comm=paste0(path_in,"/communes.shp")
comm=st_read(f_comm,stringsAsFactors = F)

# Classe de l'objet 
class(comm)
# Structure de l'objet comm
str(comm)

# Afficher la couche
plot(comm$geometry)

# Noms des champs contenus dans le df
names(comm)

# Premiers enregistrements
head(comm)

# Bbox (bounding box = emprise spatiale)
st_bbox(comm)

# Sauvegarder sous forme de shapefile
f_out=paste0(path_out,"/communes.shp") 
st_write(comm,f_out,delete_layer=TRUE) # delete_layer=T : overwrite existing file

# Sauvegarder au format geopackage
f_out=paste0(path_out,"/communes.gpkg")
st_write(comm,dsn=f_out,layer="communes",delete_layer=TRUE)

# 1 fichier geopackage peut contenir plusieurs couches
gembloux=comm[comm$ADMUNAFR=="Gembloux",] # select the polygon for Gembloux city
gembloux
st_write(gembloux,dsn=f_out,layer="gembloux",delete_layer=TRUE) # add a second layer in the gpkg file

# Afficher la liste des couches contenues dans 1 geopackage
st_layers(f_out)

# Lecture d'une couche dans 1 gpkg
f_in=paste0(path_out,"/communes.gpkg")
gembloux = st_read(dsn = f_in, layer="gembloux",stringsAsFactors = F)


# Lire les données vectorielles en format "spatvecor" (terra)
comm=terra::vect(f_in,layer="gembloux")
class(comm)

# Conversion spatvector <-> sf
comm=st_as_sf(comm) # spatvector -> sf
comm=vect(comm) # sf -> spatvector

# 3. Système de coordonnées (CRS) ------------

# Identifier le CRS
st_crs(comm)

# Définir le CRS quand il est manquant ou incorrect
f_loc=paste0(path_in,"/localites.shp")
loc=read_sf(f_loc,stringsAsFactors = F)
st_crs(loc)
# attribuer le CRS epsg:31370 à la couche "loc"
st_crs(loc)=31370
st_crs(loc)[1] # on affiche uniquement le code EPSG

# reprojeter une couche vectorielle
loc_wgs84 = st_transform(loc, 4326)
st_crs(loc_wgs84)[1]

# comparer les emprises spatiales des 2 couches
st_bbox(loc)
st_bbox(loc_wgs84)

# 4. Lecture/écriture de tables (format .csv) --------------

f_stat=paste0(path_in,"/stat_population_2018.csv")
stat_pop =read.table(f_stat,header=T,sep=";",stringsAsFactors = F)

names(stat_pop)
head(stat_pop, n = 10)
summary(stat_pop)

f_stat=paste0(path_in,"/stat_agricole_2018.csv")
stat_agri=read.table(f_stat, sep=";",dec=",", stringsAsFactors = F)
head(stat_agri)

# Sauvegarder un df dans un fichier csv
f_out=paste0(path_out,"/stat_agricole_2018.csv")
write.table(stat_agri,file=f_out,col.names = T,
            row.names = F,sep=";",dec=".")

# 5. Sélectionner/renommer des champs ---------------

f_comm=paste0(path_in,"/communes.shp")
comm=st_read(f_comm,stringsAsFactors = F)
names(comm)

# exemple 1 : conserver les 5 premiers champs de la couche "comm"
comm1 = dplyr::select(comm,OBJECTID:ADMUNAFR)
names(comm1) # le champs "geometry" est conservé par défaut

# exemple 2 : sélectionner certains champs
comm1 = dplyr::select(comm,ADMUKEY,ADMUNAFR,SHAPE_Area)
names(comm1)

# exemple 3 : sélectionner certains champs et renommer ceux-ci
comm1 = dplyr::select(comm,INS = ADMUKEY, numprov= ADPRKEY, nom = ADMUNAFR)
names(comm1)

# 6. Jointures de tables ----------------------

names(comm1) # créé au § 5
names(stat_pop) # créé au § 4

# Les champs utilisés dans la jointure doivent être du même type
class(comm1$INS)==class(stat_pop$INS)

comm2 = left_join(comm1,stat_pop, by = c("INS" = "INS"))
names(comm2)

head(comm2)
summary(comm2)

# 7. Selection par attributs ----------------

names(comm2)
plot(comm2[,"numprov"])

# Sélectionner les communes de la province de Namur
comm_nam = comm2[comm2$numprov == 9,] 
plot(comm_nam$geometry)

# 8. Selection par localisation (jointures spatiales) --------------

# Sélectionner les localités de la province de Namur (inclusion)
prov_nam=summarize(comm_nam)
loc_nam = loc[st_intersects(loc,prov_nam,sparse=F),]
plot(comm_nam$geometry)
plot(loc_nam$geometry,add=T)

# Sélectionner les localités qui sont hors de la province de Namur
loc_out_nam = loc[st_disjoint(loc,prov_nam,sparse=F),]
plot(loc_out_nam$geometry)

# Sélectionner les localités qui sont à moins de 20 km de la Centrale nucléaire de Tihange
tihange=st_read(paste0(path_in,"/tihange.gpkg"))
loc_prox_tihange = loc[st_is_within_distance(loc,tihange,20000,sparse=F),]
plot(loc$geometry)
plot(loc_prox_thiange,add=T,col="red")
plot(tihange,add=T,col="blue")

# Sélectionner les communes qui sont traversées par la rivière "Vesdre"
vesdre=st_read(paste0(path_in,"/vesdre.gpkg"))
vesdre=summarize(vesdre)

comm_vesdre = comm[st_crosses(comm,vesdre,sparse=F),]
plot(comm_vesdre$geometry)
plot(vesdre$geom,add=T,col="blue")

# 9. Création de nouveaux attributs dans 1 couche vectorielle --------------------

# 9.1 Attributs géométriques : calcul de la surface ------------------

# Surface en m²
comm_nam$surf=st_area(comm_nam)
head(comm_nam$surf)

class(comm_nam$surf)

# Convertir le résultat en km² avec 1 format "numeric"
comm_nam$surf_km2 = as.numeric(st_area(comm_nam)/1000000)
head(comm_nam$surf_km2)
class(comm_nam$surf_km2)

# 9.2 Ajouter les coordonnées des points dans le df -------------
f_loc=paste0(path_in,"/localites.shp")
loc=read_sf(f_loc,stringsAsFactors = F)
st_crs(loc)=31370
loc$x=st_coordinates(loc)[,1]
loc$y=st_coordinates(loc)[,2]
head(loc)

# Ajouter les coordonnées géographiques
loc_wgs84 = st_transform(loc,4326)
loc$long=st_coordinates(loc_wgs84)[,1]
loc$lat=st_coordinates(loc_wgs84)[,2]
head(loc)

# 9.3 Exemple de création d'un nouvel attribut -------------
#     Compléter la couche "comm" avec le nom de l'arrondissement
#     dans lequel se trouve chaque commune
#     Les données sont insérées avec 1 jointure  

# Lire le fichier arrondissement_2018.csv
f_arrond=paste0(path_in,"/arrondissements_2018.csv")
arrond=read.csv(f_arrond, sep=";" ,stringsAsFactors = F)
head(arrond)
comm_nam$INS

class(comm_nam$INS)

# step 1 : créer un champ comm_nam$id_arrond compatible avec arrond$INS
comm_nam$id_arr=as.integer(as.numeric(comm_nam$INS)/1000)*1000
comm_nam$id_arr

# step 2 : créer 1 jointure pour insérer le nom de l'arrondissement dans comm_nam
comm_nam=left_join(comm_nam,arrond,by = c("id_arr" = "INS"))

# vérifier le résultat
names(comm_nam)
plot(comm_nam[,"Arrondissement"])

# sauvegarder le résultat dans 1 gpkg
f_comm_nam=paste0(path_out,"/comm_nam.gpkg")
st_write(comm_nam,f_comm_nam,delete_layer = T)

# 10. Tableaux ou objets de synthèse avec group_by() et summarize() --------
f_comm_nam=paste0(path_out,"/comm_nam.gpkg")
comm_nam=st_read(f_comm_nam,stringsAsFactors = F)

names(comm_nam)

# calculer la densité de population par commune
comm_nam$pop_dens=comm_nam$Total/comm_nam$surf_km2
summary(comm_nam$pop_dens)

# aggréger les données par arrondissement
t1=Sys.time()
arr_nam = comm_nam %>% 
            group_by(Arrondissement) %>% 
            summarize(nb_comm = n(),
                      pop_tot = sum(Total),
                      pop_dens_moy = mean(pop_dens))
t2=Sys.time()
(t2-t1)
class(arr_nam) 
names(arr_nam)
plot(arr_nam[1]) # les polygones ont été agrégés par arrondissement

# convertir l'objet en dataframe et supprimer la géométrie des objets
# (abandon des classes "sf", "tbl_df" et "tbl")
df_comm_nam = as.data.frame(st_drop_geometry(comm_nam))
class(df_comm_nam)

t1=Sys.time()
df_arr_nam = df_comm_nam %>% 
  group_by(Arrondissement) %>% 
  summarize(nb_comm = n(),
            pop_tot = sum(Total),
            pop_dens_moy = mean(pop_dens))
t2=Sys.time()
(t2-t1)
df_arr_nam

# agréger toutes les communes de la province de Namur
prov_nam = comm_nam %>% 
  summarize(nb_comm = n(),
            pop_tot = sum(Total))
names(prov_nam)
prov_nam
plot(prov_nam$geom)

# summarize sans calcul statistique
prov_nam = comm_nam %>% 
  summarize()
names(prov_nam)

prov_nam$id=1
names(prov_nam)

# 11. Regrouper plusieurs objets sf en un seul (merge layers) ---------------

# Fusionner les shapefiles 100kmE39N29.shp et 100kmE39N30.shp
f1=paste0(path_in,"/100kmE39N30.shp")
f2=paste0(path_in,"/100kmE39N29.shp")
clc1=st_read(f1)
clc2=st_read(f2)
nrow(clc1)
nrow(clc2)

clc=rbind(clc1,clc2)
nrow(clc)

# 12. Vérifier la validité géométrique d'un objet sf --------

test=st_is_valid(clc,reason=T)
table(test)

# Corriger les erreurs de géométrie si nécessaire
clc_valid=st_make_valid(clc)
test=st_is_valid(clc_valid,reason=T)
table(test)

# 13. Intersection ---------------

# 13.1 Intersection de 2 couches (lignes et polygones) -------------

f_routes=paste0(path_in,"/routes.shp")
routes=st_read(f_routes,stringsAsFactors = F)

int=st_intersection(routes,comm_nam)
names(int)
plot(arr_nam$geom) # arr_nam : créé au § 9.3
plot(int[14],lwd=4,add=T)

# Calculer la longueur des segments de route
int$long_km = as.numeric(st_length(int))/1000 

f_out=paste0(path0,"/int_routes_comm.gpkg")
st_write(int,f_out,delete_layer=TRUE)

# Calculer la longueur de route par commune

df_int=as.data.frame(st_drop_geometry(int))
tableau1 = df_int %>% 
  group_by(nom) %>% 
  summarize(long_route = sum(long_km))
tableau1

# 13.2 Intersection de 2 couches (polygones x polygones) -----------
# Intersection des communes de la province de Namur et de la couche CLC
# Utiliser le résultat de l'intersection pour produire un tableau de surface 
# commune x classe CLC niveau 01 (CODE_01)
# CODE01 : 1=urbain, 2=agri, 3=foret, 4=z. humides, 5=eau

# Intersection de 2 couches de polygones
int=st_intersection(clc,comm_nam)
int$surf_ha=as.numeric(st_area(int))/10000

# Tableau de synthèse communes x code01
int=as.data.frame(st_drop_geometry(int)) 
head(int)

# l'abandon de la géométrie permet de réduire fortement le temps de calcul 
# pour la création du tableau 2 ci-dessous
tableau2= int %>% group_by(nom,CODE_01) %>% 
          summarize(surf_ha=sum(surf_ha))
tableau2

# Mise en forme du tableau 2 
tableau3=spread(tableau2,CODE_01,surf_ha)
names(tableau3)=c("commune","Z.urbaines","Z.agricoles","Forêt","Eau")
tableau3

# 14. Buffers ------------

# Exemple 1 : construire un buffer de 600 m autour 
# des points de la couche eoliennes.shp
eol=st_read(paste0(path_in,"/eoliennes.shp"))
eol_buff=st_buffer(eol,dist=600)
plot(eol_buff$geometry)

# Fusionner les buffers par parcs éoliens (champ [id_parc])
names(eol_buff)
parc_buff = eol_buff %>% 
  group_by(id_parc) %>% 
  summarize(nb_mat = n())
plot(parc_buff$geometry)

# Buffer de 1 km autour de la commune de Gembloux
names(comm_nam) # créé au § 7
gembloux=comm_nam[comm_nam$nom=="Gembloux",] 
buff_1km=st_buffer(gembloux,dist=1000)
plot(buff_1km$geom)
plot(gembloux$geom,add=T,col="red")

# buffer avec distance négative
buff_1km_inside=st_buffer(gembloux,dist=-1000)
plot(buff_1km_inside$geom,add=T,col="blue")

# 15. Difference  --------------------
anneau=st_difference(gembloux,buff_1km_inside)
plot(anneau$geom,col="red")

# 16. Calculs de distance ---------------

# 16.1 Recherche du plus proche voisin et calcul de distance ------------
# rechercher la ligne électrique la + proche de chaque éolienne
eol = st_read(paste0(path_in,"/eoliennes.shp"))
pow_l = st_read(paste0(path_in,"/power_lines.shp"))

# recherche de la ligne électrique la plus proche
eol$id_pow_l=st_join(eol, pow_l, join = st_nearest_feature)$full_id
eol$id_pow_l

# calcul de la distance à la ligne la plus proche
# au préalable, il faut agréger les lignes en 1 seul objet
pow_l=summarize(pow_l)
eol$d2pow_l=as.numeric(st_distance(eol,pow_l))
head(eol)

# vérifier le résultat dans QGIS
f_out=paste0(path_out,"eol_d2powl.gpkg")
st_write(eol,f_out,delete_layer=TRUE)

# 16.2 matrice de distance ------------
# chercher les 2 localités les plus proches de chaque localité
f_loc=paste0(path_in,"/localites.shp")
loc=read_sf(f_loc,stringsAsFactors = F)
st_crs(loc)=31370
list1=nngeo::st_nn(loc,loc,k=3,returnDist = T)
names(list1)
head(list1$nn)
head(list1$dist)

# Utiliser une boucle pour mettre le résultat en forme

loc$vois1=NA
loc$vois2=NA
loc$dvois1=NA
loc$dvois2=NA
# boucle pour extraire l'info
for (i in 1:nrow(loc)){
  nn=list1$nn[[i]]
  dist=list1$dist[[i]]
  loc$vois1[i]=loc$NOM[nn[2]]
  loc$vois2[i]=loc$NOM[nn[3]]
  loc$dvois1[i]=dist[2]
  loc$dvois2[i]=dist[3]
}
head(loc)
f_out=paste0(path_out,"/loc_voisins.gpkg")
st_write(loc,f_out,delete_layer = T)


# 17. Conversion entre types de géométries -----------------

# 17.1 Création de centroïdes ------------------
centro=st_centroid(comm_nam)
plot(comm_nam$geom)
plot(centro$geom,add=TRUE)

# pour imposer la présence du centroïdes à l'intérieur de l'objet
centro2=st_point_on_surface(comm_nam)
plot(centro2$geom,add=T,col="blue")

# 17.2 Convertir des polygones en lignes ou en points -----

# convertir les polygones de la couche comm en lignes
# la fonction réalise 1 conversion 1 polygone -> 1 ligne !!!!
# la couche comm étant constituée de MULTIPOLYGONS, 
# il faut au préalable transformer ceux-ci en POLYGONS simples

class(comm$geometry)
comm2=st_cast(comm,"POLYGON")
comm3=st_cast(comm2,"LINESTRING")

# vérifier dans QGIS
# la frontière séparant 2 communes est "dédoublée"
f_out=paste0(path_out,"/comm_line.gpkg")
st_write(comm3,f_out, delete_layer=TRUE)

# Convertir des lignes/polygones -> points
gembloux=comm3[comm3$ADMUNAFR=="Gembloux",]
pnt=st_cast(gembloux,"POINT")
nrow(pnt)
plot(pnt$geometry)

# éclater des lignes en segments 
segm=nngeo::st_segments(gembloux)
nrow(segm)
names(segm)
class(segm$result)
plot(segm$result)

# 17.3 Renommer le champ contenant les géométries d'1 objet sf  -----
names(segm)
names(segm)[14]="geometry"
st_geometry(segm)="geometry"
plot(segm$geometry)

# 18 Conversions  dataframe -> sf ou sf -> dataframe ---------------

# 18.1 Conversion df -> sf

# Convertir le fichier localite_wallonie.csv en 1 couche de points
f_in=paste0(path_in,"/localites_wallonie.csv")
loc_wall=read.table(f_in,header=T,sep=";",stringsAsFactors = F)
names(loc_wall)
loc_wall=st_as_sf(loc_wall,coords=c("X","Y"),crs=31370)
plot(loc_wall$geometry)
names(loc_wall)

# 18.2 Conversion sf -> df

# Supprimer la colonne de géométrie contenue dans l'objet loc_wall 
# créé au point précédent
df_loc_wall=st_drop_geometry(loc_wall)
names(df_loc_wall)

# les coordonnées ont disparu !! -> créer des champs XY au préalable
loc_wall$X=st_coordinates(loc_wall)[,1]
loc_wall$Y=st_coordinates(loc_wall)[,2]
df_loc_wall=st_drop_geometry(loc_wall)
names(df_loc_wall)

# 19 Création de géométries de base ---------------


# 19.1. générer 50 quadrats de 1 km²  répartis aléatoirement en Wallonie ----

# créer 1 couche "Wallonie"
f_in=paste0(path_in,"/communes.shp")
comm=st_read(f_in,stringsAsFactors = F)
wall=summarize(comm) # créer 1 couche "Wallonie"
plot(wall$geometry)

# step 1 : générer les centres des quadrats
centers=st_sample(comm, size=50, type = "random", exact = TRUE)  
class(centers)
centers=st_as_sf(centers)
names(centers)
names(centers)[1]="geometry"
st_geometry(centers)="geometry"
class(centers)
centers$id=seq.int(nrow(centers))
f_out=paste0(path_out,"/centers.shp")
st_write(centers,f_out,delete_layer = T)

# step 2 : générer les quadrats et les stocker dans 1 objet sf
quadrat=NULL # objet sf dans lequel vont être stockés les quadrats
for (i in 1:nrow(centers)){
  # centre du quadrat
  x0=st_coordinates(centers)[i,1]
  y0=st_coordinates(centers)[i,2]
  # création d'un objet sfg (1 géomtrie)
  sfg = st_linestring(rbind(c(x0-500, y0-500), 
                            c(x0-500, y0+500),
                            c(x0+500, y0+500),
                            c(x0+500, y0-500),
                            c(x0-500, y0-500)))
  # création d'un objet sfc (1 géomtrie + 1 CRS)
  sfc=st_sfc(sfg,crs=31370)
  # création d'un df
  df <- data.frame(id = i)
  # création d'un objet sf (sfc + df)
  sf=st_sf(cbind(sfc, df))
  # empiler les quadrats dans l'objet quadrat
  quadrat=rbind(quadrat,sf)
}

# step 3 : convertir la géométrie des objets (lignes -> polygones)
class(quadrat$geometry)
quadrat=st_cast(quadrat,"POLYGON")
plot(wall$geometry)
plot(quadrat$geometry,add=T,col="red")
class(quadrat$geometry,col="red")

# sauvegarder
f_out = paste0(path_out,"/quadrats.gpkg")
st_write(quadrat,f_out,delete_layer = T)


# 19.2. générer des ellipses représentant la zone d'influence des éoliennes ----
#       du parc d'Ernage (Gembloux) et délimiter l'emprise théorique du parc

library(DescTools) # utilisée pour appliquer une rotation aux ellipses

# Lecture des localisations des éoliennes
eol=st_read(paste0(path_in,"/eoliennes_ernage.shp"))

# paramètres de base
d_rotor = 80 # diamètre du rotor en m
a=7*d_rotor
b=4*d_rotor
azimut=235 # direction des vents dominants

emprise_eol=NULL 
# boucle sur les éoliennes
for (i in 1:nrow(eol)){

  # générer une ellipse centrée sur l'éolienne
  ell0=st_ellipse(eol[i,],a,b,res=100)
  plot(ell0)
  class(ell0)

  # orienter l'ellipse par rapport aux vents dominants
  ell0=DescTools::Rotate(st_coordinates(ell0),
                         theta=azimut/180*pi)
  class(ell0) # les données se présentent sous le forme d'une liste
  names(ell0)
  # convertir la liste -> df  
  df=as.data.frame(ell0$x)
  df
  names(df)="x"
  df$y=ell0$y
  # convertir df -> sf
  ell0=st_as_sf(df,coords=c("x","y"),crs=3812)
  plot(ell0$geometry)
  ell0$id=1

  # regrouper les points pour en faire 1 polygone
  pol = ell0 %>% 
    group_by(id) %>%
    summarize(do_union = FALSE) %>%
    st_cast("POLYGON")
  
  # empiler les résultats dans 1 sf
  emprise_eol=rbind(emprise_eol,pol)
}

# afficher le résultat
plot(emprise_eol$geometry)
plot(eol$geometry,add=T)

# tracer 1 enveloppe convexe autour des ellipses
limite = summarize(emprise_eol)
limite= st_convex_hull(limite)
plot(limite$geometry,add=T,lwd=3)

# 20. Divers -----------

# 20.1 Comblement de trous ----------------
f_in=paste0(path_in,"/compartiment.gpkg")
comp=st_read(f_in,stringsAsFactors = F)
plot(comp$geom)

comp=nngeo::st_remove_holes(comp)
plot(comp$geom)

# 20.2 polygones de voronoï-------------
# découper la Wallonie sur base de la proximité aux différentes localités
f_loc=paste0(path_in,"/localites.shp")
loc=read_sf(f_loc,stringsAsFactors = F)
st_crs(loc)=31370
# créer 1 couche avec les limites de la Wallonie
f_comm=paste0(path_in,"/communes.shp")
comm=st_read(f_comm,stringsAsFactors = F)
wall=summarize(comm) # utilisé pour découper le résulat
# créer la couche de polygones de Voronoï
loc=st_intersection(loc,wall) # élimine les localités hors Wallonie
all_loc=dplyr::summarize(loc) # regrouper les points avant de générer les pol. de voronoï
v = st_voronoi(all_loc) # création des polygones de voronoi
v = st_collection_extract(v,type="POLYGON") # convertir les géométries -> POLYGONES
v$id_vor=seq.int(nrow(v)) # attribuer 1 identifiant
# attribuer les attributs des localités aux polygones de voronoï
int=st_intersection(loc,v) # croiser les polygones avec les localités
df=as.data.frame(st_drop_geometry(int))
v=left_join(v,df, by=c("id_vor"="id_vor")) # ajouter les attributs des localités
# découper les polygones aux limites de la Wallonie
v=st_intersection(v,wall) 
plot(v$geometry)
plot(loc$geometry,add=T,col="red")
# sauvegarder les résultat
f_out=paste0(path_out,"/vor.gpkg")
st_write(v,f_out,delete_layer = T)

