1 Exercice 1
2 Exercice 2
3 Exercice 3
4 Exercice 4
5 Exercice 5
6 Exercice 6
7 Exercice 7
8 Exercice 8
9 Exercice 9
10 Exercice 10
11 Exercice 11
12 Exercice 12
13 Exercice 13
14 Exercice 14
15 Exercice 15

1 Exercice 1

Solution : on créée d’abord une fonction qui convertit les dégrés minutes secondes en dégrés décimaux. On va exprimer ces angles en radians pour pouvoir les utiliser directement avec les fonctions trigonométriques de R :

convert_degree <- function(vec) {
  # vec est un vecteur de taille 3 : c(deg, minute, second)
  deg <- vec[1]
  minute <- vec[2] 
  second <- vec[3]
  stopifnot(minute < 60 & second < 60)
  minute <- minute + second / 60
  (deg + minute / 60) * pi / 180
}

Ensuite, on créée une fonction qui calcule la distance entre deux points :

dist_deg <- function(coord_A, coord_B, R = 6378) {
  # convert
  long_A <- convert_degree(coord_A[[1]]) 
  lat_A <- convert_degree(coord_A[[2]]) 
  long_B <- convert_degree(coord_B[[1]]) 
  lat_B <- convert_degree(coord_B[[2]]) 
  # distance
  R * acos(sin(lat_A) * sin(lat_B) + cos(lat_A) * cos(lat_B) * cos(long_A - long_B))
}

Application :

dist_deg(list(c(2, 21, 7), c(48, 51, 24)),
         list(c(1, 26, 3), c(43, 36, 16)))

## [1] 588.9053

2 Exercice 2

En utilisant le package photon ou tidygeocoder:

my_address <- c("compans cafarelli, toulouse", "INRAE, Auzeville")
locgeo.full <- photon::geocode(my_address, limit = 1, key = "place")
tidygeocoder::geo(address = my_address, method = "osm", 
                   lat = latitude, long = longitude)

## # A tibble: 2 x 3
##   address                     latitude longitude
##   <chr>                          <dbl>     <dbl>
## 1 compans cafarelli, toulouse     43.6      1.44
## 2 INRAE, Auzeville                43.5      1.50

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

3 Exercice 3

On montre ici le code utilisé.

Pour la France, on a utilisé le CRS officiel (EPSG:2154).
Pour la carte du monde, on a utilisé la projection de World Robinson. Pour d’autres représentations pour la surface globale de la terre, on pourra consulter cette page : https://mgimond.github.io/Spatial/coordinate-systems-in-r.html ainsi que des représentations originales sur : http://magrit.cnrs.fr/docs/projections_fr.html
Pour des représentations de l’Europe, cette page : https://observablehq.com/@toja/five-map-projections-for-europe

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)

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)

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)

4 Exercice 4

Pour obtenir l’abcisse : \(10000+5\times 500 + \frac{1}{2}\times500=12750\) et pour l’ordonnée : \(15000+6\times 500=18000\)

5 Exercice 5

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")

En utilisant la norme sf:

library(sf)
city_sf <- st_as_sf(city, 
                 coords = c("Longitude","Latitude"), 
                 crs = 4326) # On indique le systeme de coordonnées utilisé, ici WGS84

En utilisant la norme sp:

library(sp)
city_sp <- city
coordinates(city_sp) <- ~Longitude+Latitude
proj4string(city_sp) <- CRS(SRS_string = "EPSG:4326")

library(leaflet)
m <- leaflet(city_sf) %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addCircles()
m

6 Exercice 6

library(sf)
bike <- read_sf(dsn = "Donnees/velo toulouse/velo-toulouse.shp")
iris <- read_sf(dsn = "Donnees/velo toulouse/pop.shp")

st_crs(bike)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## 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["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

st_crs(iris)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## 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["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

7 Exercice 7

library(raster)
pop_fr <- raster("Donnees/pop_fr/fra_pd_2020_1km.tif")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

8 Exercice 8

On utilise la syntaxe à la ggplot

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.0     ✓ dplyr   1.0.5
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()

ggplot(iris, aes(geometry = geometry)) +
  geom_sf()

On utilise la syntaxe classique :

plot(st_geometry(iris))

9 Exercice 9

iris_tlse <- iris %>%
  group_by(code_insee) %>%
  summarise(p12_pop = sum(p12_pop))

## although coordinates are longitude/latitude, st_union assumes that they are planar

plot(st_geometry(iris_tlse))

iris_tlse

## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1.350379 ymin: 43.53266 xmax: 1.515341 ymax: 43.66866
## Geodetic CRS:  WGS 84
## # A tibble: 1 x 3
##   code_insee p12_pop                                                    geometry
##   <chr>        <dbl>                                               <POLYGON [°]>
## 1 31555      453317. ((1.366204 43.59024, 1.36626 43.59038, 1.365661 43.5908, 1…

10 Exercice 10

plot(st_geometry(iris_tlse), lwd = 2)
plot(st_geometry(iris), add = T, lty = 2)
plot(st_geometry(bike), add = T, pch = 16, cex = 0.5, col = "royalblue")

11 Exercice 11

Les jeux de données iris et bike sont exprimés dans le WGS 84.

st_crs(bike)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## 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["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

st_crs(iris)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## 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["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

On va les transformer dans le système officiel en France dont le EPSG vaut 2154

bike <- st_transform(bike, crs = 2154)
iris <- st_transform(iris, crs = 2154)
ggplot(data = iris, aes(geometry = geometry)) +
   geom_sf() + 
   geom_sf(data = bike) +
   ggtitle("Velo Toulouse")

plot(st_geometry(iris), axes = T)
plot(st_geometry(bike), add = T, pch = 16, col = "royalblue")

12 Exercice 12

iris$area <- st_area(iris)

13 Exercice 13

Dans la fonction st_intersection(), on positionne d’abord l’objet bike et ensuite l’objet iris :

bike <- st_intersection(bike, iris)

14 Exercice 14

On trouve d’abord la localisation géographique du domicile et du lieu de travail :

address <- c("2 impasse Armand Carrel, 31000 Toulouse, France",
             "1 Esplanade de l'université, 31000 Toulouse")
locgeo.full <- photon::geocode(address, limit = 1, key = "place")
locgeo.full

##                                          location     osm_id osm_type name
## 1 2 impasse Armand Carrel, 31000 Toulouse, France 2025525717        N <NA>
## 2     1 Esplanade de l'université, 31000 Toulouse 7216063191        N <NA>
##   housenumber                    street postcode     city     state country
## 1           2     Impasse Armand Carrel    31200 Toulouse Occitanie  France
## 2           1 Esplanade de l'Université    31000 Toulouse Occitanie  France
##   osm_key osm_value      lon      lat  msg
## 1   place     house 1.465337 43.63364 <NA>
## 2   place     house 1.435015 43.60464 <NA>

Ensuite, on calcule le chemin optimal :

library(osrm)

## Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright

## Routing: OSRM - http://project-osrm.org/

short_path <- osrmRoute(
    src = locgeo.full[1, c("lon", "lat")],
    dst = locgeo.full[2, c("lon", "lat")],
    returnclass = "sp",
    overview = "full"
  )

On importe le fonds de carte :

xy <- SpatialPoints(locgeo.full[, c("lon", "lat")], 
                    CRS("EPSG:4326")) 
library(cartography)

## Warning in fun(libname, pkgname): rgeos: versions of GEOS runtime 3.9.0-CAPI-1.16.2
## and GEOS at installation 3.8.0-CAPI-1.13.1differ

mtqOSM <- getTiles(x = xy, type = "osm", zoom = 13)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

On représente le trajet :

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("My home", "My Work"))

15 Exercice 15

La réponse a été inspirée de ce post : https://rcarto.github.io/Intro_a_la_carto_avec_R/ex/ville.html

On convertit les données sur les villes dans le CRS de Robinson :

city_sf <- st_transform(city_sf, "ESRI:54030")

On représente la carte avec les villes :

plot(st_geometry(world_az))
plot(st_geometry(city_sf), add = T, pch = 6, cex = 0.2)

par(mar = c(0,0,1.2,0), bg = "white")
plot(st_geometry(world_az), col = "lightblue", border = "white", lwd = 0.2, bg = "cornsilk2")
# ajout d'un ombrage pour les pays du monde
plot(st_geometry(world_az) + c(100000,-100000), add=T, col = "grey30", border = NA)
plot(st_geometry(world_az), col = "ivory4", border ="ivory3", lwd = 0.5, add=TRUE)
# choix des bornes de la discrétisation du taux de croissance
bks <- c(min(city_sf$X2015.2020), 0, 1, 2, 3, 4, max(city_sf$X2015.2020))
# choix des couleurs
cols <- carto.pal("wine.pal", 1, "green.pal", 5)
propSymbolsChoroLayer(x = city_sf, var = "X2020", border = "white", legend.title.cex = 0.6,
                 lwd = 0.5, symbols = "circle", inches = 0.2,
                 var2 = "X2015.2020", breaks = bks, col = cols,
                 legend.var.title.txt = "Population en 2020\n(en milliers)",
                 legend.var2.title.txt = "Taux de croissance\nannuel moyen\n(2015-2020)",
                 legend.var.pos = "topright", 
                 legend.var2.pos = "topleft", legend.values.cex = 0.6)
layoutLayer(title = "Population des agglomérations urbaines en 2018", 
            sources = "WUP, 2018", author = "Inspiré de T. Giraud, 2017", col = "black", 
            coltitle = "white", scale = NULL, frame = FALSE,  tabtitle = TRUE)

Correction exercice

Thibault Laurent

Mise à jour : 03 mai 2021