## Cartographie avec R
# Formation INRAE, 4 Mai 2021
# Thibault.Laurent@tse-fr.eu

## ---- eval = F------------------------------------------------------------------------------
## install.packages(c(
##   "cartography", # réaliser des cartes
##   "classInt", # discrétisation de variables quantitatives
##   "geogrid", # permet de transformer un polygone en hexagone
##   "ggspatial", # syntaxe complémentaires à la ggplot
##   "GISTools", # outils pour faire de la carto
##   "leaflet", # interractivité avec JavaScript
##   "mapview", # interractivité avec JavaScript
##   "maptools", # manipulation de données "spatial",
##   "OpenStreetMap", # OSM
##   "osrm", # openstreetmap avec R
##   "raster", # manipulation de données raster
##   "RColorBrewer", # palette de couleurs pour carto
##   "readxls", # import de données XLS
##   "rgdal", # import de données spatiales
##   "rgeos", # manipulation de données spatiales
##   "sf", # nouvelle classe d'objets spatials
##   "sp", # ancienne classe d'objets spatials
##   "tidygeocoder", # localiser une adresse
##   "tidyverse", # ggplot, dplyr, etc
##   "tmaptools" # pour la carto
##   ),
##                  dependencies = TRUE)
## devtools::install_github(repo = 'rCarto/photon') # calcul de distance routière
## devtools::install_github("rCarto/popcircle")


## ---- eval = F------------------------------------------------------------------------------
## folder <- getwd()
## download.file("http://www.thibault.laurent.free.fr/cours/spatial/Donnees.zip",
##   destfile = paste0(folder, "/Donnees.zip"))
## unzip("Donnees.zip", exdir = paste0(folder, "/Donnees"))


## ----sig, echo=FALSE, fig.cap="Exemple de SIG.", out.width = '70%', fig.align = "center"----
knitr::include_graphics("./Figures/columbus2.jpg")


## ----carto_issue, echo=FALSE, fig.cap="Deux fichiers, deux représentations.", out.width = '100%', fig.align = "center"----
knitr::include_graphics("./Figures/referentiel.jpeg")


## +proj=longlat +datum=WGS84 +no_defs


## "EPSG:4326"


## GEOGCRS["WGS 84",

##     DATUM["World Geodetic System 1984",

##         ELLIPSOID["WGS 84",6378137,298.257223563,

##             LENGTHUNIT["metre",1]]],

##     PRIMEM["Greenwich",0,

##         ANGLEUNIT["degree",0.0174532925199433]],

##     CS[ellipsoidal,2],

##         AXIS["geodetic latitude (Lat)",north,

##             ORDER[1],

##             ANGLEUNIT["degree",0.0174532925199433]],

##         AXIS["geodetic longitude (Lon)",east,

##             ORDER[2],

##             ANGLEUNIT["degree",0.0174532925199433]],

##     USAGE[

##         SCOPE["Horizontal component of 3D system."],

##         AREA["World."],

##         BBOX[-90,-180,90,180]],

##     ID["EPSG",4326]]


## ----geoide, echo=FALSE, fig.cap="Le géoïde ou la véritable forme de la terre.", out.width = '40%', fig.align = "center"----
knitr::include_graphics("./Figures/geoide.jpg")


## ----ellipsoid, echo=FALSE, fig.cap="Ajustement d'un ellipsoïde au géoïde.", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/ellipsoid_shape.jpg")


## -------------------------------------------------------------------------------------------
rgdal::projInfo("ellps")


## ----Clarke, echo=FALSE, fig.cap="Les datum Clarke 1880 et WGS 84.", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/i10_EllipsoideClarkeGRS80Certu.jpg")


## ----local, echo=FALSE, fig.cap="Exemple de datum local.", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/ellipsoid_local2.JPG")


## DATUM["World Geodetic System 1984",

##         ELLIPSOID["WGS 84",6378137,298.257223563,

##             LENGTHUNIT["metre",1]]]


## +datum=WGS84 +ellps=WGS84


## DATUM["Nouvelle_Triangulation_Francaise_Paris",

##             SPHEROID["Clarke 1880 (IGN)",6378249.2,293.4660212936265,

##                 AUTHORITY["EPSG","7011"]],

##             TOWGS84[-168,-60,320,0,0,0,0],

##                 AUTHORITY["EPSG","6807"]]


## +a=6378249.2 +b=6356515.000000472 +towgs84=-168,-60,320,-0,-0,-0,0


## ----coord, echo=FALSE, fig.cap="Exemple de datum local.", out.width = '30%', fig.align = "center"----
knitr::include_graphics("./Figures/coordgeo.JPG")


## -------------------------------------------------------------------------------------------
require("photon")
address <- c("21 allee de Brienne, 31000 Toulouse, France",
             "2 Rue du Doyen Gabriel Marty, 31042 Toulouse")
locgeo.full <- photon::geocode(address, limit = 1, key = "place")
locgeo.full


## -------------------------------------------------------------------------------------------
library(leaflet)
m <- leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addMarkers(lng = locgeo.full$lon, lat = locgeo.full$lat, popup = "MT")
m


## -------------------------------------------------------------------------------------------
osm_s1 <- tidygeocoder::geo(
     address = address, method = "osm",
     lat = latitude, long = longitude
)


## ----world, echo=FALSE, fig.cap="Exemple de représentation globale", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/world.jpeg")


## ----cyl, echo=FALSE, fig.cap="Projection cylindrique (vue globale)", out.width = '30%', fig.align = "center"----
knitr::include_graphics(c("./Figures/Projection_cylindrique.jpg",
                          "./Figures/Projection_cylindrique2.jpg"))


## ---- eval = T------------------------------------------------------------------------------
sf::st_crs("+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs")


## ----utm31, echo=FALSE, fig.cap="Projection UTM en France utilisée sur les cartes de l'IGN", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/utm_31.jpeg")


## ----conic, echo=FALSE, fig.cap="Projection conique (vue globale et tronquée)", out.width = '40%', fig.align = "center"----
knitr::include_graphics("./Figures/conic.jpg")
knitr::include_graphics("./Figures/conic2.jpg")


## ---- eval = T------------------------------------------------------------------------------
sf::st_crs("+proj=aea +lat_1=29.83 +lat_2=45.83 +lat_0=37.5 
     +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")


## ----albers, echo=FALSE, fig.cap="Projection Albers", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/albers.jpeg")


## ---- eval = T------------------------------------------------------------------------------
sf::st_crs("+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 
     +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs")


## ----lambert, echo=FALSE, fig.cap="Projection Lambert 93", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/lambert.jpeg")


## ----azimuth, echo=FALSE, fig.cap="Projection azimuthale", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/azimuth.jpg")


## ----azimuth2, echo=FALSE, fig.cap="Exemple de projection azimuthale : US National Atlas Equal Area", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/azimuthal.jpeg")


## ---- eval = T------------------------------------------------------------------------------
sf::st_crs("EPSG:2154")


## -------------------------------------------------------------------------------------------
EPSG <- rgdal::make_EPSG()


## -------------------------------------------------------------------------------------------
EPSG[grep("RGF93", EPSG$note), 1:2]


## ---- echo = F, message = F, warning = F, fig.width=10, fig.height = 6----------------------
library(sf)
world_sf <- read_sf("Donnees/World WGS84/Pays_WGS84.shp") 
fr <- world_sf[53, ]
fr_2 <- st_transform(fr, 2154)
op <- par(mfrow=c(1,2))
plot(st_geometry(fr), axes = T, main = "EPSG:4326")
plot(st_geometry(fr_2), axes = T, main = "EPSG:2154")
par(op)
title("Proposition 1", line = 3)


## ---- echo = F, message = F, warning = F, fig.width=10, fig.height = 12---------------------
world_az <- st_transform(world_sf, st_crs('ESRI:54030'))
op <- par(mfrow=c(2,1))
plot(st_geometry(world_sf), axes = T, main = "EPSG:4326")
plot(st_geometry(world_az), axes = T, main = "ESRI:54030")
par(op)
title("Proposition 2", line = 3)


## ---- echo = F, message = F, warning = F, fig.width=10, fig.height = 6----------------------
europe <- world_sf[c(8, 9,10,11,12,13,14, 15,16,17,18,19,20,22,24,25,26,
                     27, 28, 31, 32, 33, 34, 35, 36, 37, 40, 41,43, 46, 47, 48, 49, 50,
                     52,53, 57, 56, 58, 59, 61, 64, 71, 74, 80, 85, 93),]
europe_az <- st_transform(europe, st_crs('ESRI:102031'))
op <- par(mfrow=c(1, 2))
plot(st_geometry(europe), axes = T, main = "EPSG:4326")
plot(st_geometry(europe_az), axes = T, main = "ESRI:102031")
par(op)
title("Proposition 3", line = 3)


## ----vec_rast, echo=FALSE, fig.cap="Vecteur VS raster", out.width = '70%', fig.align = "center"----
knitr::include_graphics("./Figures/vector_vs_raster.jpg")


## ----raster, echo=FALSE, fig.cap="Exemple de fichier raster", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/spatial4.jpeg")


## ---- echo = F------------------------------------------------------------------------------
plot(0, 0, xaxt = "n", yaxt = "n", xlim = c(0, 10), ylim = c(0, 10),
     xlab = "", ylab = "")
text(0,0, "O", pos = 2)
points(5.5, 6, pch = 16, col = "red")
abline(v=0:10)
abline(h=0:10)


## ----vecteur, echo=FALSE, fig.cap="Exemple de fichier vecteur", out.width = '100%', fig.align = "center"----
knitr::include_graphics("./Figures/250px-Simple_vector_map.jpg")


## {"type": "FeatureCollection",

##   "crs": {"type": "name", "properties": {"name": "EPSG:26904"}},

##   "features": [

##     {"type": "Feature",

##       "properties": {

##         "name": "Van Dorn Street",

##         "marker-color": "#0000ff",

##         "marker-symbol": "rail-metro",

##         "line": "blue"},

##       "geometry": {

##         "type": "Point",

##         "coordinates": [

##           -77.12911152370515,

##           38.79930767201779

##         ]}

##     }]

## }


## ---- eval = F------------------------------------------------------------------------------
## download.file(url = "https://raw.githubusercontent.com/riatelab/basemaps/master/World/countries.geojson",
##               destfile = "Donnees/country.geojson")
## download.file(url = "https://raw.githubusercontent.com/riatelab/basemaps/master/World/graticule30.geojson",
##               destfile = "Donnees/graticule.geojson")


## -------------------------------------------------------------------------------------------
seisme_df <- read.csv2("Donnees/earthquake/earthquakes.csv")
head(seisme_df, 2)


## -------------------------------------------------------------------------------------------
library(sp)
seisme_sp <- seisme_df
coordinates(seisme_sp) <- ~Longitude + Latitude 
class(seisme_sp)


## -------------------------------------------------------------------------------------------
str(seisme_sp)


## -------------------------------------------------------------------------------------------
head(seisme_sp@data, 2)
head(seisme_sp@coords, 2)


## -------------------------------------------------------------------------------------------
proj4string(seisme_sp) <- CRS(SRS_string = "EPSG:4326")


## ---- eval = F------------------------------------------------------------------------------
## proj4string(seisme_sp) <- CRS("+proj=longlat +datum=WGS84 +no_defs
##                             +ellps=WGS84 +towgs84=0,0,0")
## # équivalent à cette ligne de commande
## proj4string(seisme_sp) <- CRS("+proj=longlat +datum=WGS84 +no_defs")


## -------------------------------------------------------------------------------------------
library(sf)
seisme_sf <- st_as_sf(seisme_df, coords = c("Longitude", "Latitude"))
class(seisme_sf)


## -------------------------------------------------------------------------------------------
head(seisme_sf)


## ----sf class, echo=FALSE, fig.cap="Simple feature geometry", out.width = '50%', fig.align = "center"----
knitr::include_graphics("./Figures/sf-classes.png")


## -------------------------------------------------------------------------------------------
head(st_drop_geometry(seisme_sf), 2)


## -------------------------------------------------------------------------------------------
seisme_sf %>% 
  st_drop_geometry() %>%
  head(2)


## -------------------------------------------------------------------------------------------
st_crs(seisme_sf) <- 4326


## ---- eval = T------------------------------------------------------------------------------
library(leaflet)
m <- leaflet(seisme_sf) %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addCircles()
m


## ---- eval = F------------------------------------------------------------------------------
## # Population data base
## download.file(url = "https://population.un.org/wup/Download/Files/WUP2018-F12-Cities_Over_300K.xls",
##               destfile = "Donnees/population_mondiale/citypop.xls")
## # Growth rate data base
## download.file(url = "https://population.un.org/wup/Download/Files/WUP2018-F14-Growth_Rate_Cities.xls",
##               destfile = "Donnees/population_mondiale/cityevo.xls")


## ---- message = F---------------------------------------------------------------------------
library(readxl)
citypop <- data.frame(read_excel("Donnees/population_mondiale/citypop.xls", skip = 16))
cityevo <- data.frame(read_excel("Donnees/population_mondiale/cityevo.xls", skip = 16))
city <- merge(x = citypop, y = cityevo[, c(4,9:24)], by = "City.Code")


## -------------------------------------------------------------------------------------------
library(rgdal)
world_sp <- readOGR(dsn = "Donnees/World WGS84", 
                    layer = "Pays_WGS84")


## -------------------------------------------------------------------------------------------
str(world_sp[53, ])


## -------------------------------------------------------------------------------------------
str(world_sp[53, ]@polygons[[1]])


## ---- warning = F---------------------------------------------------------------------------
cat(wkt(world_sp))


## -------------------------------------------------------------------------------------------
world_sf <- read_sf("Donnees/World WGS84/Pays_WGS84.shp") 
st_crs(world_sf)


## -------------------------------------------------------------------------------------------
head(world_sf)


## -------------------------------------------------------------------------------------------
str(st_geometry(world_sf[53, ]))


## -------------------------------------------------------------------------------------------
library(raster)
wind <- raster("Donnees/wind/FRA_wind-speed_200m.tif")


## -------------------------------------------------------------------------------------------
wind


## -------------------------------------------------------------------------------------------
wind[2000:2001, 4000:4002, 1]
wind[2000:2001, 4000:4002]


## ---- eval = F------------------------------------------------------------------------------
## leaflet() %>% addTiles() %>%
##   addRasterImage(wind)


## -------------------------------------------------------------------------------------------
methods(class = "Spatial")
methods(class = "sf")
methods(class = "raster")


## -------------------------------------------------------------------------------------------
dim(world_sp)
dim(world_sf)


## -------------------------------------------------------------------------------------------
row.names(world_sp) <- as.character(world_sp@data$NOM)


## -------------------------------------------------------------------------------------------
maghreb_sp <- world_sp[c("Algeria", "Morocco",
                      "Libya", "Mauritania", 
                      "Tunisia", "Western Sahara"),]


## ---- message = F---------------------------------------------------------------------------
library(tidyverse)
maghreb_sf <- world_sf %>%
  filter(NOM %in% c("Algeria", "Morocco", "Libya",
 "Mauritania", "Tunisia", "Western Sahara"))


## -------------------------------------------------------------------------------------------
plot(maghreb_sp, col = RColorBrewer::brewer.pal(6, "Set3"))
title("Maghreb countries")
legend("topleft", legend = row.names(maghreb_sp), cex = 0.6,
 title = "Nom des pays", fill = RColorBrewer::brewer.pal(6, "Set3"))


## -------------------------------------------------------------------------------------------
plot(maghreb_sf, main = "Maghreb countries")


## ---- warning = F---------------------------------------------------------------------------
library(maptools)
maghreb_sp_new <- unionSpatialPolygons(maghreb_sp, 
              IDs = c("Algeria", "Morocco", "Libya", 
                      "Mauritania", "Tunisia", "Morocco"))
class(maghreb_sp_new)


## -------------------------------------------------------------------------------------------
maghreb2.df <- data.frame(NOM = c("Algeria", "Morocco", "Libya",
                                 "Mauritania", "Tunisia"), 
                          pib = c(4177, 2875, 4337, 1139, 3861), 
                          pop = c(39728, 34663, 6418, 4046, 11179),
                          region = rep("N", 5))
row.names(maghreb2.df) <- maghreb2.df$NOM


## -------------------------------------------------------------------------------------------
maghreb_sp_new <- SpatialPolygonsDataFrame(maghreb_sp_new, maghreb2.df) 


## -------------------------------------------------------------------------------------------
sapply(maghreb_sp_new@polygons, function(x) x@ID)


## -------------------------------------------------------------------------------------------
maghreb_sp_new$pib_tot <- maghreb_sp_new$pib * maghreb_sp_new$pop
big_maghreb <- aggregate(maghreb_sp_new[, c("pib_tot", "pop")], 
                         by = list(name = rep("maghreb", 5)), 
          FUN = sum)
big_maghreb$pib <- round(big_maghreb$pib_tot / big_maghreb$pop)
maghreb_sp_new$pib_tot <- NULL


## ---- message = F---------------------------------------------------------------------------
maghreb_sf_new <- maghreb_sf %>%
  group_by(NOM = c("Tunisia", "Algeria", "Morocco", 
                   "Libya", "Morocco", "Mauritania")) %>%
  summarise() %>%
  merge(maghreb2.df, by = "NOM")


## -------------------------------------------------------------------------------------------
egypt_sp <- world_sp["Egypt", ]
egypt_df <- data.frame(NOM = "Egypt", pib = 3599, pop = 92442, region = "N")
egypt_sp <- merge(egypt_sp, egypt_df, by = "NOM")
row.names(egypt_sp) = "Egypt"


## ---- message = F, warning = F--------------------------------------------------------------
northAf_sp <- spRbind(maghreb_sp_new, egypt_sp)


## -------------------------------------------------------------------------------------------
egypt_sf <- world_sf[world_sf$NOM == "Egypt", ]
egypt_sf <- merge(egypt_sf, egypt_df)
northAf_sf <- maghreb_sf_new %>%
  rbind(egypt_sf)


## -------------------------------------------------------------------------------------------
plot(world_sp)
plot(seisme_sp[1:1000, ], pch = 16, cex = 1, add = TRUE)
title("Strongest earthquake in 40 years")


## -------------------------------------------------------------------------------------------
plot(st_geometry(world_sf))
plot(st_geometry(seisme_sf)[1:1000, ], pch = 16, cex = 1, add = TRUE)
title("Strongest earthquake in 40 years")


## -------------------------------------------------------------------------------------------
ggplot(data = world_sf, aes(geometry = geometry)) +
   geom_sf() + 
   geom_sf(data = seisme_sf[1:1000,]) +
   ggtitle("Strongest earthquake in 40 years")  


## -------------------------------------------------------------------------------------------
plot(world_sp, col = rgb(0.95, 0.95, 0.95), bg="#A6CAE0", 
     xlim = c(-20, 40), ylim = c(15, 40), lty = 2)
plot(maghreb_sp_new, col = rgb(1, 1, 0.3), add = T)
title("GDP in North Africa in 2012")
text(coordinates(maghreb_sp_new),
 paste(row.names(maghreb_sp_new), maghreb_sp_new$pib, sep = " : "), cex = 0.8)
legend("bottomright","Source : CIA", cex = 0.6)


## ---- warning = F---------------------------------------------------------------------------
ggplot(data = world_sf, aes(geometry = geometry)) +  
  geom_sf() + 
  geom_sf(data = northAf_sf, fill = rgb(1, 1, 0.3)) +
  geom_sf_text(data = northAf_sf, 
     aes(label = paste(NOM, pib, sep = " : ")), 
      size = 3) +
  coord_sf(xlim = c(-20, 40), ylim = c(15, 40)) +
  ggtitle("GDP in North Africa in 2012", 
    subtitle = "Source : CIA") 


## -------------------------------------------------------------------------------------------
france_sp <- world_sp["France", ]
my_crs_fr <- CRS(SRS_string = "EPSG:2154")
france_sp2 <- spTransform(france_sp, my_crs_fr)


## -------------------------------------------------------------------------------------------
my_crs_fr
cat(wkt(my_crs_fr))


## -------------------------------------------------------------------------------------------
france_sf <- world_sf %>%
  filter(NOM == "France")
france_sf2 <- st_transform(france_sf, my_crs_fr)
france_sf2 <- st_transform(france_sf, 2154)


## -------------------------------------------------------------------------------------------
st_crs(france_sf2)


## -------------------------------------------------------------------------------------------
crs_robi <- CRS("ESRI:54030")


## -------------------------------------------------------------------------------------------
world_sp_robi <- spTransform(x = world_sp, CRSobj = crs_robi)
seisme_sp_robi <- spTransform(seisme_sp, CRSobj = crs_robi)


## -------------------------------------------------------------------------------------------
world_sf_robi <- st_transform(x = world_sf, crs = crs_robi)
seisme_sf_robi <- st_transform(seisme_sf, crs = crs_robi)
ggplot(data = world_sf_robi, aes(geometry = geometry)) +
   geom_sf() + 
   geom_sf(data = seisme_sf_robi[1:1000,]) +
   ggtitle("Strongest earthquake in 40 years")  


## -------------------------------------------------------------------------------------------
b0 <- st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))


## -------------------------------------------------------------------------------------------
b1 <- b0 + 2
b2 <- b0 + c(-0.2, 2)
x <- st_sfc(b1, b2)
plot(b0, axes = T, xlim = c(-1, 3), ylim = c(-2, 4))
plot(x, border = 'red', axes = T, add = T)


## -------------------------------------------------------------------------------------------
a0 <- b0 * 0.8
a1 <- b0 * 1.2 
a2 <- b0 * 0.8 + c(2, 0.7)
y <- st_sfc(a0,a1,a2)
plot(b0, axes = T, xlim = c(-1, 3), ylim = c(-2, 4))
plot(y, border = 'green', add = T)


## -------------------------------------------------------------------------------------------
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
c1 <- b0 * rot(pi/4)
c2 <- b0 * rot(pi/3) + c(2, -0.5)
c3 <- (b0 * 2) * rot(pi/3) + c(2, -0.5)
z <- st_sfc(c1, c2, c3)
plot(b0, axes = T, xlim = c(-1, 3), ylim = c(-2, 4))
plot(z, border = 'magenta', add = T)


## ---- message = F---------------------------------------------------------------------------
library(rgeos)
world_sp_robi$area <- gArea(world_sp_robi, byid = T)
world_sf_robi$area <- st_area(world_sf_robi)


## ---- echo = F------------------------------------------------------------------------------
plot(c(0.25, 0.5, 1.5), c(0.25, 0.5, 0.5), xlim = c(0, 2), ylim = c(0, 1),
     xaxt = "n", yaxt = "n", pch = 16, col = "grey", xlab = "", ylab = "")
text(c(0.25, 0.5, 1.5), c(0.25, 0.5, 0.5), c("s1", "s2", "s3"), pos = 2)
text(c(0.2, 1.2), c(1, 1), c("T1", "T2"), col = "blue", pos = 1)
segments(0,0,2,0, col = "blue")
segments(0,1,2,1, col = "blue")
segments(0,0,0,1, col = "blue")
segments(1,0,1,1, col = "blue")
segments(2,0,2,1, col = "blue")


## ---- warning = F---------------------------------------------------------------------------
japan_sp <- world_sp["Japan",]


## -------------------------------------------------------------------------------------------
bb.jp_sp <- as(extent(bbox(japan_sp)), "SpatialPolygons")


## -------------------------------------------------------------------------------------------
proj4string(bb.jp_sp) <- CRS(SRS_string = "EPSG:4326")


## -------------------------------------------------------------------------------------------
ind <- over(seisme_sp, bb.jp_sp)


## -------------------------------------------------------------------------------------------
seisme_sp.jp <- seisme_sp[which(!is.na(ind)), ]


## -------------------------------------------------------------------------------------------
plot(japan_sp, border = NA, col = "NA", bg="#A6CAE0", axes = T)
plot(world_sp, col = "#E3DEBF", border = NA, add = T)
plot(seisme_sp.jp, add = TRUE, col = "red", cex = 0.7)


## -------------------------------------------------------------------------------------------
japan_sf <- world_sf %>% filter(NOM == "Japan")


## -------------------------------------------------------------------------------------------
bb.jp_sf <- st_as_sfc(st_bbox(japan_sf))


## -------------------------------------------------------------------------------------------
seisme.jp_sf <- st_intersection(seisme_sf, bb.jp_sf) 


## -------------------------------------------------------------------------------------------
seisme.jp_sf[, "Magnitude"] %>% ggplot(aes(geometry = geometry)) +  
  geom_sf(data = world_sf) + 
  geom_sf(aes(colour = Magnitude)) +
  xlim(125, 145) + 
  ylim(30, 48) + 
  theme_bw()


## -------------------------------------------------------------------------------------------
crs_jp <- CRS("EPSG:2443")
seisme_sp.jp2 <- spTransform(seisme_sp.jp, crs_jp)
japan_sp2 <- spTransform(japan_sp, crs_jp)


## -------------------------------------------------------------------------------------------
japan_sp_buff <- gBuffer(japan_sp2, width = 50000)


## -------------------------------------------------------------------------------------------
seisme_sp.jp3 <- gIntersection(seisme_sp.jp2, japan_sp_buff)


## -------------------------------------------------------------------------------------------
plot(japan_sp2, border = NA, col = "NA", bg="#A6CAE0", axes = T)
plot(japan_sp2, col = "#E3DEBF", border = NA, add = T)
plot(japan_sp_buff, add = TRUE, lty = 2)
plot(seisme_sp.jp3, add=T, col = "red", cex = 0.7)


## -------------------------------------------------------------------------------------------
seisme.jp_sf2 <- st_transform(seisme.jp_sf, 2443)
japan_sf2 <- st_transform(japan_sf, 2443)


## -------------------------------------------------------------------------------------------
japan_sf_buff <- st_buffer(japan_sf2, dist = 50000)


## ---- warning = F---------------------------------------------------------------------------
seisme.jp_sf3 <- st_intersection(seisme.jp_sf2, japan_sf_buff)


## -------------------------------------------------------------------------------------------
ggplot(data = world_sf, aes(geometry = geometry)) +
  geom_sf(bg = "#A6CAE0", col = "#E3DEBF") +
  geom_sf(data = japan_sf_buff, lty=2, col = scales::alpha("green", 0.1)) +
  geom_sf(data = japan_sf2, fill = "#A6CAE0") + 
  geom_sf(data = seisme.jp_sf3, col = "red", cex = 0.7) +
  xlim(125, 145) + ylim(30, 48) + theme_bw()


## ---- message = F---------------------------------------------------------------------------
suneu <- raster("Donnees/soleil/pvgis_g13year00.asc", 
 crs = CRS(SRS_string = "EPSG:4326"), native = TRUE)


## -------------------------------------------------------------------------------------------
suneu[150, 685, 1]


## -------------------------------------------------------------------------------------------
plot(suneu)


## -------------------------------------------------------------------------------------------
suneu.bb <- crop(suneu, extent(bbox(france_sp))) 


## -------------------------------------------------------------------------------------------
suneu.fr <- rasterize(france_sp, suneu.bb, mask = TRUE)


## -------------------------------------------------------------------------------------------
plot(france_sp)
image(suneu.fr, add = TRUE)


## -------------------------------------------------------------------------------------------
dim(suneu)


## -------------------------------------------------------------------------------------------
suneu.agg <- aggregate(suneu, fact = c(15, 10), fun = mean)


## -------------------------------------------------------------------------------------------
dim(suneu.agg)


## ---- message = F, warning = F--------------------------------------------------------------
library(GISTools)
plot(france_sp)
image(suneu.agg, add = TRUE)
masker <- poly.outer(suneu.agg, france_sp)
add.masking(masker)


## ---- eval = T------------------------------------------------------------------------------
france_dep <- getData("GADM", country = "france", level = 2)  


## ---- eval = T, warning = F-----------------------------------------------------------------
france_alt <- getData("alt", country = "FRA", mask = TRUE)


## ---- echo = F, eval = F--------------------------------------------------------------------
## save(france_dep, france_alt, file = "./Donnees/france_dep.RData")


## ---- echo = F, eval = F--------------------------------------------------------------------
## load("./Donnees/france_dep.RData")


## -------------------------------------------------------------------------------------------
plot(france_alt)
plot(france_dep, add = T)


## ---- message = F---------------------------------------------------------------------------
long <- c(1.432022, 1.436879, 1.435409)
lat <- c(43.605417, 43.606688, 43.610372)
xy <- SpatialPoints(cbind(long, lat), 
                    CRS("EPSG:4326")) 


## ---- warning = F---------------------------------------------------------------------------
library(cartography)
mtqOSM <- getTiles(x = xy, type = "osm", zoom = 13)


## -------------------------------------------------------------------------------------------
st_crs(mtqOSM)


## -------------------------------------------------------------------------------------------
tilesLayer(mtqOSM)
title("UT1-C")
plot(xy, col = 2, cex = 2, add = TRUE)
text(coordinates(xy), col = 2, pos = 3,
  labels = c("Manu", "Arsenal","Metro"))   


## ---- eval = T------------------------------------------------------------------------------
library(osrm)
short_path <- osrmRoute(
    src = xy[1, ],
    dst = xy[2, ],
    returnclass = "sp",
    overview = "full"
  )


## ---- echo = F, eval = F--------------------------------------------------------------------
## save(short_path, file = "short_path.RData")


## ---- echo = F, eval = F--------------------------------------------------------------------
## load("short_path.RData")


## -------------------------------------------------------------------------------------------
tilesLayer(mtqOSM)
title("Short path")
plot(xy, col = 2, cex = 2, add = TRUE)
plot(short_path, add = T, col = "royalblue", lwd = 2)
text(coordinates(xy), col = 2, pos = 3,
  labels = c("Manu", "Arsenal","Metro"))   


## ----bertin1, echo=FALSE, fig.cap="Extrait 1 de l'ouvrage de Bertin", out.width = '100%', fig.align = "center"----
knitr::include_graphics("./Figures/bertin_1.jpg")


## ----bertin2, echo=FALSE, fig.cap="Extrait 2 de l'ouvrage de Bertin", out.width = '100%', fig.align = "center"----
knitr::include_graphics("./Figures/bertin_2.jpg")


## -------------------------------------------------------------------------------------------
pays.af <- c("Seychelles", "Equatorial Guinea", "Gabon", "Botswana", "Mauritius", 
 "South Africa", "Namibia", "Angola", "Swaziland", "Congo", "Cape Verde", "Ghana", 
 "Sudan", "Djibouti", "Nigeria", "Sao Tome and Principe", "Cameroon", "Lesotho", 
 "Gambia, The", "Chad", "Senegal", "Kenya", "Ivory Coast", "Zambia", "Burkina Faso",
 "Tanzania, United Republic of", "Benin", "Rwanda", "Uganda", "Comoros", "Guinea-Bissau", 
 "Mali", "Mozambique", "Guinea", "Ethiopia", "Madagascar", "Malawi", "Togo", 
 "Sierra Leone", "Niger", "Central African Republic", "Eritrea", "Burundi", "Somalia", 
 "Zimbabwe", "Liberia", "Zaire")


## -------------------------------------------------------------------------------------------
africa.df <- data.frame(pib = c(14745, 11279, 7382, 6800, 9260, 5735, 5033, 4167, 3689, 
 1761, 3043, 1766, 2487, 2676, 2730, 1596, 1327, 1112, 661, 776, 1219, 1337, 1426, 
 1338, 575, 947, 784, 728, 709, 1242, 603, 751, 590, 769, 641, 467, 381, 571, 588, 
 361, 377, 715, 306, 293, 1445, 710, 497),
                     pop = c(93, 1169, 1948, 2121, 1263, 55386, 2315, 27884, 1104, 
 4856, 525, 27849, 38903, 914, 181137, 199, 23298, 2059, 2086, 14111, 14578, 47878, 
 23226, 15879, 18111, 51483, 10576, 11369, 38225, 777, 1737, 17439, 27042, 11432, 
 100835, 24234, 16745, 7323, 7172, 20002, 4493, 4475, 10160, 13797, 13815, 4472, 76245), 
                     region = c("E", "C", "C", "Au", "E", "Au", "Au", "C", "Au", "C",
 "O", "O", "E", "E", "O", "C", "C", "Au", "O", "C", "O", "E", "O", "E", "O", "E", "O",
 "E", "E", "E", "O", "O", "E", "O", "E", "E", "E", "O", "O", "O", "C", "E", "E", "E",
 "E", "O", "C"), NOM = pays.af)	


## -------------------------------------------------------------------------------------------
africa.sub_sp <- merge(world_sp[pays.af, ], africa.df, by = "NOM")
row.names(africa.sub_sp) <- as.character(africa.sub_sp$NOM)
africa.sub_sf <- merge(world_sf %>% filter(NOM %in% pays.af), africa.df, by = "NOM") 


## ---- warning = F---------------------------------------------------------------------------
africa_sp <- spRbind(northAf_sp, africa.sub_sp)
africa_sf <- rbind(northAf_sf, africa.sub_sf)


## ---- fig.width=12, fig.height = 8----------------------------------------------------------
require(RColorBrewer)
display.brewer.all(type="seq")


## ---- fig.width=12, fig.height = 8----------------------------------------------------------
display.brewer.all(type="qual") 


## ---- fig.width=12, fig.height = 8----------------------------------------------------------
display.brewer.all(type="div")


## -------------------------------------------------------------------------------------------
pal.reg <- RColorBrewer::brewer.pal(5, "Set3") 


## -------------------------------------------------------------------------------------------
ind_reg <- as.numeric(factor(africa_sp$region))


## -------------------------------------------------------------------------------------------
crs_af <- CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m")
africa_sp <- spTransform(africa_sp, crs_af)


## -------------------------------------------------------------------------------------------
scalebar <- function(loc, length, unit = "km", division.cex = 0.8, 
                     bg = "white", border = "black",...) {

  x <- c(0, length / 2, length) + loc[1]
  y <- c(0, length / 10, length / 9.5, length / 5.5) + loc[2]
  cols <- c("black", "white")
  for (i in 1:2) rect(x[i], y[1], x[i+1], y[2], col = cols[i])
  for (i in 1:3) segments(x[i], y[2], x[i], y[3])
  labels <- round(c(0, length / 2, length) / 1000)
  labels[3] <- paste0(labels[3], " km") 
  text(x[c(1, 2, 3)], y[4], labels, adj = 0.5, cex = division.cex)
}

north.arrow <- function(x, y, h, lab = "North", lab.pos = "below") {
  polygon(c(x, x, x + h/2), c(y - (1.5*h), y, y - (1 + sqrt(3)/2) * h), 
          col = "black", border = NA)
  polygon(c(x, x + h/2, x, x - h/2), c(y - (1.5*h), y - (1 + sqrt(3)/2) * h, 
                                       y, y - (1 + sqrt(3)/2) * h))
  if(lab.pos == "below") text(x, y - (2.5 * h), lab, adj = c(0.5, 0), cex = 1)
  else text(x, y + (0.25 * h), lab, adj = c(0.5, 0), cex = 1.5)
}


## -------------------------------------------------------------------------------------------
plot(africa_sp, col = pal.reg[ind_reg])
title("Africa's regions")
legend("topright", legend = c("Aust", "Cent", "East", "North", "West"),
       cex = 0.8, title = "Region", fill = pal.reg)
north.arrow(6*10^6, -2*10^6, 5*10^5, lab = "N")
scalebar(c(-2*10^6, -4*10^6), 2500000)


## -------------------------------------------------------------------------------------------
library(ggspatial)
ggplot(data = africa_sf, aes(geometry = geometry)) +
  geom_sf(aes(fill = as.factor(region)), legend="meh") +
  scale_fill_discrete("Regions") +
  coord_sf(crs = st_crs(crs_af)) +
    annotation_scale(location = "bl", width_hint = 0.5) +
    annotation_north_arrow(location = "bl", which_north = "true", 
        pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
        style = north_arrow_fancy_orienteering)


## -------------------------------------------------------------------------------------------
ggplot(data = africa_sf, aes(geometry = geometry)) +
  geom_sf(aes(fill = pib)) +
  coord_sf(crs = st_crs(crs_af)) +
    annotation_scale(location = "bl", width_hint = 0.5) +
    annotation_north_arrow(location = "bl", which_north = "true", 
        pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
        style = north_arrow_fancy_orienteering)


## -------------------------------------------------------------------------------------------
pib <- africa_sp$pib
library(classInt)
library(RColorBrewer)
pal1 <- brewer.pal(5, "Blues")
opar <- par(mfrow = c(2, 2), mar = c(3, 2, 1, 0), las = 1)
plot(classIntervals(pib, 5, "equal"), pal = pal1, main = "equal")
plot(classIntervals(pib, 5, "quantile"), pal = pal1, main = "quantile")
plot(classIntervals(pib, 5, "kmeans"), pal = pal1, main = "K-means")
plot(classIntervals(pib, 5, "jenks"), pal = pal1, main = "Jenks")
par(opar)


## -------------------------------------------------------------------------------------------
bk <- round(classIntervals(pib, 5, "kmeans")$brks, digits = 1)
ind <- findInterval(pib, bk, all.inside = TRUE)
plot(africa_sp, col = pal1[ind])
decoup <- c("<=1678.5", "]1678.5;3321]", "]3321;6267.5]", 
            "]6267.5;10269.5]", ">10269.5")
legend("bottomleft", legend = decoup, cex = 0.8, title = "GDP", 
       fill = pal1)
SpatialPolygonsRescale(layout.north.arrow(1), offset = c(50, -30), 
                       scale = 5, plot.grid = F)
title("GDP per inhabitant")


## ---- message = F, warning = F--------------------------------------------------------------
library(geogrid)
africa_cells_hex <- calculate_grid(shape = africa_sf, grid_type = "hexagonal", 
                                   seed = 3, verbose = F)
africa_hex <- assign_polygons(africa_sf, africa_cells_hex) %>%
  st_set_crs(2154)


## -------------------------------------------------------------------------------------------
plot(st_geometry(africa_hex), col = pal1[ind])
decoup <- c("<=1678.5", "]1678.5;3321]", "]3321;6267.5]", 
            "]6267.5;10269.5]", ">10269.5")
legend("bottomleft", legend = decoup, cex = 0.8, title = "GDP", 
       fill = pal1)
SpatialPolygonsRescale(layout.north.arrow(1), offset = c(50, -30), 
                       scale = 5, plot.grid = F)
title("GDP per inhabitant")


## -------------------------------------------------------------------------------------------
choroLayer(spdf = africa_sp, var = "pib", method = "fisher-jenks", legend.pos = "topleft")
layoutLayer(title = "GDP in Africa", sources = "World Bank", frame = TRUE, 
            col = "NA",  scale = NULL)


## ---- warning = F---------------------------------------------------------------------------
plot(africa_sp, border = NA, col = "NA", bg = "#A6CAE0")
plot(world_sp, col = "#E3DEBF", border = NA, add = T)
plot(africa_sp, col = "#D1914D", border = "grey80", add = T)
propSymbolsLayer(spdf = africa_sp, var = "pib" )


## ---- warning = F---------------------------------------------------------------------------
library(popcircle)
res <- popcircle(x = africa_sf, var = "pop")
circles <- res$circles
shapes <- res$shapes
par(mar = c(0, 0, 0, 0))
plot(st_geometry(circles), col = "#bcd39c", 
     border = "white", bg = "#eafdcf")
plot(st_geometry(shapes), col = "#fffc99", 
     border = "#fffc99", add = TRUE)
labelLayer(x = circles[1:20,], txt = "NOM",
           halo = TRUE, col ="#8e8358", 
           cex = seq(.8,.4, length.out = 20),
           font = 2, bg = "white", r = .15, 
           overlap = FALSE)


## ---- eval = F------------------------------------------------------------------------------
## library(leaflet)
## m <- leaflet() %>%
##   addTiles() %>%  # Add default OpenStreetMap map tiles
##   addMarkers(lng = 1.432022, lat = 43.605417, popup = "MT")


## ---- eval = F------------------------------------------------------------------------------
## library(mapview)
## mapview(seisme.jp_sf)


## ---- message = F---------------------------------------------------------------------------
library("spatstat")
library("maptools")
af <- aggregate(africa_sp, list(fr = rep("Af", nrow(africa_sp))), FUN = length)
S <- as(af, "owin")


## -------------------------------------------------------------------------------------------
af_ppp <- as.ppp(coordinates(africa_sp), W = S)


## -------------------------------------------------------------------------------------------
marks(af_ppp) <- africa_sp$pib[-25]


## -------------------------------------------------------------------------------------------
op <- par(mfrow = c(2, 2), mar = c(0, 0, 3.3, 0))
plot(af_ppp, main = "PIB")
plot(nnmark(af_ppp), main = "Plus proche voisin")
plot(Smooth(af_ppp), main = "Lissage spatial")
plot(idw(af_ppp, power = 2), main = "Lissage par l'inverse de la distance")
par(op)