# Manipulation de données spatiales dans R avec le package`SF`


Le package `sf` (Spatial Feature) est dédié à la manipulation, la transformation et l'analyse de données spatiales. A la manière du Tidyverse, il combine les fonctionnalités de `sp`, `rgeos` et `rgdal` dans un package.

> Une attention particulière a été portée à la compatibilité du package avec la syntaxe pipe et les opérateurs du tidyverse.
>

install.packages("sf")

library(sf)

Ressources :

* https://r-spatial.github.io/sf/
* https://geocompr.robinlovelace.net/geometric-operations.html
* https://rcarto.github.io/carto_avec_r/chapitre1.html

***Cheatsheet*** > https://github.com/rstudio/cheatsheets/blob/master/sf.pdf


___

## Importer un shapefile

ardt <- st_read(dsn = "C:/Users/Xo/Desktop/GeoDATA/arrondissements.shp", stringsAsFactors = FALSE) 

metro <- st_read(dsn = "C:/Users/Xo/Desktop/GeoDATA/metroParis.shp", stringsAsFactors = FALSE) 

IRIS <- st_read(dsn = "C:/Users/Xo/Desktop/GeoDATA/iris.shp", stringsAsFactors = FALSE) 

velos <- st_read(dsn = "C:/Users/Xo/Desktop/GeoDATA/reseau-cyclable.shp", stringsAsFactors = FALSE) 

toilettes <- st_read(dsn = "C:/Users/Xo/Desktop/GeoDATA/sanisettesparis.shp", stringsAsFactors = FALSE) 

___

## Visualiser les arrondissements

De base toute les variables sont cartographiées

plot(ardt)


Plot en fonction d'une variable

plot(ardt["l_ar"])


## Visualiser les IRIS

Cartographier les IRIS en fonction de l'arrondissement

plot(IRIS["nom_com"])


___

## Ecrire un shapefile

st_write(ardt, "C:/Users/Xo/Desktop/GeoDATA/arrondissements2.shp")

## Ecrire un geopackage avec plusieurs couches

st_write(metro, "C:/Users/Xo/Desktop/metro.gpkg", "metro", append= TRUE)

st_write(ardt, "C:/Users/Xo/Desktop/metro.gpkg", "ardt", append= TRUE)

___
## Voir le SCR d'une couche

st_crs(ardt)


___
## Reprojeter une couche

ardt2154<- st_transform(ardt, 2154)

___
## Géoférencer un CSV

**Géoréférencer le csv des arbres**

* Importer le csv des arbres
* Créer les deux colonnes lat/long
* Tranformer le en objet SF
* Cartographier les arbres en fonction de l'arrondissement

Arbres <- Arbres %>% separate(col= geo_point_2d, 
into = c("Latitude", "Longitude"), sep= ",")

arbresgeo = st_as_sf(arbres, coords=c("Longitude","Latitude"), crs=4326)

plot(arbresgeo["ARRONDISSEMENT"])


**Géoréférencer le csv des emplacements de stationnement de vélos**


## Jointure attributaire

* Importer le fichier csv ParisINSEE
* Faire la jointure
* cartographier la variable de médiane du niveau de vie (MNV)

ArdtINSEE <- merge(ardt2154, ParisINSEE, by.x = c("c_arinsee"), by.y = c("CodeArdt"))


* Pour changer la discrétisation il suffit d'ajouter un paramètre "breaks"

* Pour changer la position de la légende il suffit d'ajouter un paramètre dans la commande (1=below, 2=left, 3=above and 4=right)

plot(ArdtINSEE["MNV"], breaks = "quantile", key.pos = 1)


## Centroides

Afficher les centroïdes des IRIS

plot(st_centroid(IRIS["iris"]))

Créer une couche des centroïdes des IRIS

 centroIRIS <- st_centroid(IRIS)

---

## Operateurs topographiques

### Calcul d'aire

Ajouter à la couche des arrondissements une colonne surface

ardtok <- ardtok %>% mutate(surface2 = st_area(ardtok))

Calculer et cartographier les surfaces des IRIS

IRIS <- IRIS %>% mutate(surface = st_area(IRIS))
plot(IRIS["surface"])

___
### Calcul de perimètre

install.packages("lwgeom")
library(lwgeom)

Ajouter à la couche arrondissement une colonne périmètre

ardtok <- ardtok %>% mutate(longueur2 = st_perimeter(ardtok))

___

### Calcul de longueur

Cartographier la couche vélos en fonction du statut des tronçons


Ajouter à la couche vélo une colonne longueur

N'oublier pas de reprojeter cette couche en 2154 !

velos <- st_transform(velos, 2154)
velos <- velos %>% mutate(longueur = st_length(velos))

Cartographier la nouvelle variable en utilisant les quantiles

plot(velos["longueur"], breaks = "quantile", key.pos = 1)


___
### Faire un **buffer**

Faire un buffer de 200m autour des stations de métro

Buffermetro <- st_buffer(metro, dist=200)

plot(Buffermetro["ligne"])


Regrouper les buffers

Buffermetro <-st_union(Buffermetro)


Faire une zone d'influence de 300m autour des toilettes publiques


___
### Dissolve

Regrouper les IRIS par arrondissement

ARDT <- IRIS %>% group_by(insee_com) %>% summarise()

___

### Dissolve with stat

Regrouper les IRIS par arrondissement avec un résumé statistique de la surface (somme)

ARDT <- IRIS %>% group_by(insee_com) %>% 
summarise(NbIris = n(), Surface = sum(surface))

---

## Jointure spatiale

**Ajouter à chaque station de métro le nom de l'arrondissement**

spatialjoin <- st_join(metro, ardt2154["l_ar"])

plot(spatialjoin["l_ar"])


**Ajouter à chaque arbres le nom de l'IRIS**


## Compter le nombre de point dans un polygone

**Compter le nombre de metro par IRIS**

IRIS <- IRIS %>% mutate(nbMetro = lengths(st_intersects(IRIS, metro)))

**Compter le nombre d'arbres par IRIS**

IRIS <- IRIS %>% mutate(nbarbres = lengths(st_intersects(IRIS, arbresgeo)))

**Compter à la fois le nombre de toilettes et de stations de métro dans chaque Arrondissement**

Ardt <- Ardt %>% mutate(nbMetro = lengths(st_intersects(Ardt, metro)),
nbtoilettes = lengths(st_intersects(Ardt, WC)))


**Compter pour chaque IRIS le nombre de stations de métros, le nombre d'arbres, le nombre de toilettes et le nombre de stationnements de vélos**


**Garder uniquemebt les IRIS contenant :**
* plus de 2 stations de métro
* plus de 500 arbres
* plus de 3 emplacements de vélos


## Jointure spatiale avec résumé statistique (somme)

**Calculer la somme totale des transactions DVF 2017 par arrondissement**

SumDVFARDT <- ardt2154 %>% st_join(DVF) %>%
group_by(l_ar) %>% summarise( total = sum(prix))

**Calculer les longueurs de piste cyclables par IRIS**

IRISvelo <- IRIS %>% st_join(velos) %>% 
     group_by(nom_iris) %>% 
     summarise( longueurvelo = sum(longueur))
     

plot(IRISvelo["longueurvelo"], breaks = "quantile", key.pos = 1)
     


## Jointure spatiale avec résumé statistique (moyenne)

**Calculer la moyenne des prix au m² par IRIS et faite une carte**

M2moyenIRIS <- IRIS %>% st_join(DVF) %>% 
group_by(nom_iris) %>% 
summarise( prixm2moyen = mean(M2OK))


**Calculer la surface moyenne des biens par IRIS et cartographier cette nouvelle variable**


___
## Calculer un ensemble de statistiques issues de l'agrégation

DVFardt <- ardt2154 %>% st_join(DVF) %>% 
group_by(l_ar) %>% 
summarise( nb = n(), total = sum(prix), prixmoyen = mean(prix), prixm2 = mean(M2OK) )

___
## Créer une grille

**Créer une grille de type hexagonale de 500m**

grille <- st_make_grid(
   ardt2154,
   cellsize = 500,
   crs = 2154, 
   what = "polygons",
   square = FALSE)


Pour changer la forme agir sur le paramètre `square` True ou False

grille2 <- st_make_grid(
   ardt2154,
   cellsize = 500,
   crs = 2154, 
   what = "polygons",
   square = TRUE)


**Compter le nombre d'arbres dans la grille hexagonale**

 grille <- st_sf(grille)
 grille <- grille %>% mutate(nbarbres = lengths(st_intersects(grille, arbresgeo)))


**Compter le nombre de transactions DVF dans la grille hexagonale**


**Calculer le prix moyen au m² dans le carroyage hexagonale**

Il faut d'abor créer un champ de regroupement comme une ID puis procéder à l'agrégation spatiale avec résumé statistique

grille <- grille %>% mutate(id = row_number())

grille <- grille %>% mutate(IDOK = id)

grillem2 <- grille %>% st_join(DVF) %>% 
    group_by(IDOK) %>%
    summarise( prixm2moyen = mean(M2OK))

plot(grillem2["prixm2moyen"], breaks = "quantile", key.pos = 1)


**Mettre en place un petit script qui :**
* Crée la grille à la volée
* Fait l'agrégation du prix moyen au m²
* Sort une carte

grille300 <- st_make_grid(
   ardt2154,
   cellsize = 300,
   crs = 2154, 
   what = "polygons",
   square = TRUE)
      
grille <- st_sf(grille300)
 
grille <- grille %>% mutate(id = row_number())

grille <- grille %>% mutate(IDOK = id)

grillem2 <- grille %>% st_join(DVF) %>% 
    group_by(IDOK) %>%
    summarise( prixm2moyen = mean(M2OK))

plot(grillem2["prixm2moyen"], breaks = "quantile", key.pos = 1)


## Intersection

**Ne garder que les stations de métro du 1er arrondissement**

Premierardt <- ardt2154 %>% filter(c_arinsee == 75101)
metro1er <- st_intersection(metro, Premierardt)

**Ne garder que les arbres du 16eme arrondissement**

Seizerardt <- ardt2154 %>% filter(c_arinsee == 75116)
arbres16eme <- st_intersection(arbresgeo, Seizerardt)


**Produire une couche avec tous les transactions DVF situés à moins de 100m d'une station de métro**

Buffermetro <- st_buffer(metro, dist=100)
DVFmetro <- st_intersection(DVF, Buffermetro)


**Garder tous les arrêts de métro situés à moins de 100m d'une piste cyclable**

pistecyclables <- velos %>% filter(typologie_s == 'Bandes cyclables')

metroprevelos <- st_intersection(metro, st_buffer(pistecyclables, 100)) %>% distinct(id_ref_zdl)

**Produire une couche avec tous les arbres situés à moins de 200m d'un métro et 100m d'un toilette publique**

Buffermetro <- st_buffer(metro, dist=200)
BufferWC <- st_buffer(WC, dist=100)
buffer <- st_intersection(Buffermetro,BufferWC)
Arbresendanger <- st_intersection(arbresgeo, buffer)


**Garder tous les toilettes situés à moins de 100m d'une piste cyclable et à moins de 100m d'un métro**

**Produire une couche avec toutes les transaction DVF situées à moins de 300m d'un métro, à plus de 100m d'un toilette publique et à moins de 100m d'une piste cyclable**

Buffermetro <- st_buffer(metro, dist=300)
Buffervelo <- st_buffer(velos, dist=100)
BufferWC <- st_buffer(WC, dist=100)
BufferWC <- st_union(BufferWC)
Inclusion <- st_intersection(Buffermetro,Buffervelo)
Inclusion <- st_union(Inclusion)
Exclusion <- st_difference(Inclusion, BufferWC)
DVFOK <- st_intersection(DVF, Exclusion)

## Trouver l'entité la plus proche

## Jouer avec les Trips

### Importer le dataset

tracesoiseaux <- st_read(dsn = "D:/Cours/13-Automne2020/DataViz_Cergy/Data_KeplerGL/9_Traces_Oiseaux.geojson")

### Transformer l'objet SF en dataframe

tracesoiseaux <- as.data.frame(tracesoiseaux)

### Changer le timestamp en temps POSIX

tracesoiseaux$timestamp <- as_datetime(tracesoiseaux$timestamp)
tracesoiseaux <- tracesoiseaux %>% mutate(UNIX= as.numeric(as.POSIXct(timestamp)))

### Créer une colone altitude (si elle n'existe pas)

tracesoiseaux<- tracesoiseaux %>% mutate(altitude = 10)

### Concaténer la colonne altitude et de temps dans une seule colonne avec un séparateur facile à transformer par la suite

tracesoiseaux <- tracesoiseaux %>% mutate(other=paste(altitude, UNIX, sep='12345'))

### Créer le Geojson (sous forme de ligne) avec la structuration attendue

 Tripsoiseaux <- tracesoiseaux %>%
     st_as_sf(coords = c("long", "lat", "other")) %>%
     group_by(tag_ident) %>%
     summarise(do_union = FALSE) %>%
     st_cast("LINESTRING")

### Ecrire le GeoJson

st_write(Tripsoiseaux, "tripoiseaux.geojson", append = TRUE)

### Remplacer dans le GeoJson le séparateur bidon (12345) par une ,