R Spatial¶

Il existe de nombreuses bibliothèques R permettant de traiter de l'information géographique. Néanmoins, je déconseille fortement d'en utiliser de nombreuses, car un petit nombre de bibliothèques est largement suffisant pour effectuer les principaux traitements géographiques. De plus, l'utilisation de nombreuses bibliothèques multiplie les dépendances et favorise l'obsolescence de vos programmes. Utiliser de nombreuses bibliothèques R révèle généralement une mauvaise maitrise de R. Ainsi, la bibliothèque sf, présentée dans le précèdent notebook pour charger des données géographiques, permet d'effectuer les principaux géotraitements proposés par les SIG, pas besoin d'aller chercher ailleurs.

Ainsi, après avoir fait un peu le tour de la bibliothèque sf, seules les bibliothèques spdep et sp seront présentées. En effet, ces dernières permettent alors d'effectuer des tâches que même les SIG classiques ne proposent pas forcément.

Géotraitement¶

Les données nécessaires pour cette partie sont disponibles ici : https://sergelhomme.fr/data/Donnees.zip

Comme vu dans le notebook précédent, la bibliothèque sf (et sa fonction st_read()) permet de lire et donc charger dans R les principaux formats de données géographiques.

In [ ]:
library(sf)
In [2]:
ad <- file.choose()
dep <- st_read(ad)
plot(dep$geometry)

Reading layer `Departements' from data source 
  `C:\Users\Serge\Downloads\Donnees(4)\Departements.shp' using driver `ESRI Shapefile'
Simple feature collection with 96 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 99227.52 ymin: 6049649 xmax: 1242375 ymax: 7110525
Projected CRS: RGF93 / Lambert-93
No description has been provided for this image
In [3]:
print(dep[1:10,])

Simple feature collection with 10 features and 16 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 522391.1 ymin: 6164368 xmax: 1077506 ymax: 7009225
Projected CRS: RGF93 / Lambert-93
   ID_GEOFLA CODE_DEPT                NOM_DEPT              NOM_CHF CODE_REG
1          1        01                     AIN      BOURG-EN-BRESSE       82
2          2        02                   AISNE                 LAON       22
3          3        03                  ALLIER              MOULINS       83
4          4        04 ALPES-DE-HAUTE-PROVENCE      DIGNE-LES-BAINS       93
5          5        05            HAUTES-ALPES                  GAP       93
6          6        06         ALPES-MARITIMES                 NICE       93
7          7        07                 ARDECHE               PRIVAS       82
8          8        08                ARDENNES CHARLEVILLE-MEZIERES       21
9          9        09                  ARIEGE                 FOIX       73
10        10        10                    AUBE               TROYES       21
                   NOM_REGION Pop_hommes Pop_femmes Pop_Etrang Pop_totale
1                 RHONE-ALPES     287336     294023      45109     581355
2                    PICARDIE     262006     276794      13447     538790
3                    AUVERGNE     164392     178416       9394     342807
4  PROVENCE-ALPES-COTE-D'AZUR      76580      81387       6388     157965
5  PROVENCE-ALPES-COTE-D'AZUR      65494      68714       4192     134205
6  PROVENCE-ALPES-COTE-D'AZUR     511716     572712     100388    1084428
7                 RHONE-ALPES     151678     159776       9343     311452
8           CHAMPAGNE-ARDENNE     139346     144856      11382     284197
9               MIDI-PYRENEES      73313      76894       6864     150201
10          CHAMPAGNE-ARDENNE     146609     154726      14327     301327
   Naissances  Deces Solde_nat Solde_mig Tx_Etrange Superfice
1       61000  37114     23886     41991     0.0776  5784.556
2       62669  48087     14582    -11105     0.0250  7419.598
3       30652  38838     -8186      6378     0.0274  7380.145
4       13741  13972      -231     18513     0.0404  6995.771
5       12882  10793      2089     10485     0.0312  5689.465
6      102510 101409      1101     71461     0.0926  4295.464
7       29910  28091      1819     23473     0.0300  5565.608
8       32094  25080      7014    -12941     0.0400  5243.878
9       12445  15900     -3455     16309     0.0457  4906.472
10      31923  25772      6151      3077     0.0475  6025.735
                         geometry
1  MULTIPOLYGON (((919195 6541...
2  MULTIPOLYGON (((735604.5 68...
3  MULTIPOLYGON (((753768.9 65...
4  MULTIPOLYGON (((992637.8 63...
5  MULTIPOLYGON (((1012912 640...
6  MULTIPOLYGON (((1018256 627...
7  MULTIPOLYGON (((831641.3 63...
8  MULTIPOLYGON (((842093.3 69...
9  MULTIPOLYGON (((631546.3 61...
10 MULTIPOLYGON (((796594.5 67...

Pour traiter de l'information géographique, il convient de maitriser des fonctions permettant de mesurer des objets géométriques pour calculer des superficies, des périmètres, des longueurs ou encore des distances. Ainsi, la fonction st_area() permet de mesurer simplement l'aire de polygones à condition de lui donner des géométries compatibles, comme celles obtenues à partir de la fonction st_read()...

In [4]:
area <- st_area(dep$geometry)
print(area)

Units: [m^2]
 [1]  5773959873  7426966023  7366255429  6994778675  5684618604  4295917034
 [7]  5560374059  5249567884  4916098055  6020256031  6356187589  8771018731
[13]  5093274202  5605610547  5768116966  5963536397  6901617185  7295951970
[19]  5891135415  8786956827  6974746097  5588515046  9211279557  5244387692
[25]  6550483841  6035864147  5927471939  6736808054  4036424863  4722538867
[31]  5876082075  6366180433  6304297299 10142248954  6235032711  6835817110
[37]  6884963482  6148759791  7865220722  5042300647  9351804263  6409620191
[43]  4797401174  4998614179  6899411553  6805778932  5223836955  5383600746
[49]  5170840751  7218483478  6004678208  8189630079  6250744280  5206508442
[55]  5280374315  6236563702  6860153216  6251518890  6862936838  5767423305
[61]  5893638150  6141503409  6714202204  8000254765  7695210885  4529303978
[67]  4153030131  4791887450  3527587589  3254724503  5381923326  8597244900
[73]  6239046366  6262545023  4595227255   105341755  6327046796  5923214834
[79]  2305306911  6030106988  6205986203  5781526792  3728154752  6037299426
[85]  3579108995  6759512956  7022591763  5550602881  5891514376  7449495200
[91]   609935151  1818648924   174911956   237807018   244891032  1252980367

Ces mesures de superficies sont associées à des unités de mesure. Si vous souhaitez récupérer uniquement les valeurs, vous pouvez utiliser la fonction as.numeric().

In [5]:
dep$Sup <- as.numeric(area) / 1000000
print(dep[c('Sup','Superfice')])

Simple feature collection with 96 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 99227.52 ymin: 6049649 xmax: 1242375 ymax: 7110525
Projected CRS: RGF93 / Lambert-93
First 10 features:
        Sup Superfice                       geometry
1  5773.960  5784.556 MULTIPOLYGON (((919195 6541...
2  7426.966  7419.598 MULTIPOLYGON (((735604.5 68...
3  7366.255  7380.145 MULTIPOLYGON (((753768.9 65...
4  6994.779  6995.771 MULTIPOLYGON (((992637.8 63...
5  5684.619  5689.465 MULTIPOLYGON (((1012912 640...
6  4295.917  4295.464 MULTIPOLYGON (((1018256 627...
7  5560.374  5565.608 MULTIPOLYGON (((831641.3 63...
8  5249.568  5243.878 MULTIPOLYGON (((842093.3 69...
9  4916.098  4906.472 MULTIPOLYGON (((631546.3 61...
10 6020.256  6025.735 MULTIPOLYGON (((796594.5 67...

La fonction st_length() permet de calculer la longueur de lignes. On pourra l'utiliser pour calculer les périmètres de polygones, à condition de récupérer les contours de ces polygones avec la fonction st_cast() par exemple.

In [6]:
st_length(dep$geometry)
st_length(st_cast(dep, "MULTILINESTRING"))

Units: [m]
 [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[39] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[77] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Units: [m]
 [1]  477704.17  597363.80  538629.71  572888.88  571079.57  412550.56
 [7]  431612.95  449016.05  487188.15  466045.65  518297.87  586191.92
[13]  565728.58  493931.53  514853.93  499165.52  772645.78  531960.72
[19]  452901.43  575904.17  651350.35  436083.55  581242.71  491534.37
[25]  639619.16  520644.04  472594.39 1044733.28  501482.80  413959.19
[31]  594503.86  776768.72  505633.16  711472.44  482349.87  556002.46
[37]  470675.55  463501.85  609568.32  467933.39  565073.41  578107.01
[43]  513300.62  457748.79  596500.39  519781.23  437443.33  473001.69
[49]  392070.74  556987.06  623667.21  568831.42  535666.24  436204.00
[55]  682508.87  506445.53  850108.04  665554.10  486418.00  717389.03
[61]  587047.14  555869.63  607035.25  507597.89  605345.29  492057.46
[67]  385279.46  461349.54  310746.68  399271.17  468609.86  653328.38
[73]  450790.79  513173.79  401574.95   50171.94  455454.84  495126.19
[79]  293748.15  554183.35  565086.51  459570.66  507083.55  535145.41
[85]  392921.18  570174.89  558455.31  491773.36  553854.47  582625.82
[91]  146522.82  249975.47   80608.49   93085.45   90969.61  247320.35

La fonction st_centroid() permet de calculer les centroides de polygones. Si l'on veut que les centroides soient situés obligatoirement à l'intérieur des polygones, on pourra utiliser la fonction st_point_on_surface(). On pourra superposer dans un graphique deux couches (les départements et leur centroide) en effectuant successivement deux plots, mais il faudra utiliser l'argument add dans le deuxième plot pour ne pas créer deux graphiques indépendants.

In [7]:
centres <- st_centroid(dep$geometry)
plot(dep$geometry)
plot(centres, add=TRUE)
No description has been provided for this image

La fonction st_distance() permet de calculer des distances minimales entre tous les types d'objets géométriques.

In [8]:
st_distance(dep$geom[1],dep$geom[2]) # Distance minimale entre les départements de l'Ain et de l'Aisne
st_distance(centres[1],centres[2]) # Distance entre les centroides de l'Ain et de l'Aisne

Units: [m]
         [,1]
[1,] 280977.5
Units: [m]
         [,1]
[1,] 407157.5

On peut utiliser les centroides pour effectuer des cartes avec des symboles proportionnels.

In [9]:
plot(dep$geometry, col="#EEEEEE")
plot(centres, add=TRUE, cex=sqrt(dep$Pop_totale) / 500, pch = 21, bg="#EEA861", lwd = 1, col = "black")
No description has been provided for this image

Pour réaliser une carte, il faudra au minimum rajouter une légende. Celle pour les symboles proportionnels est un peu complexe à réaliser lorsque l'on débute en R. Néanmoins, une fois que l'on a compris comment on peut faire, on pourra réutiliser le code pour d'autres cartes, les modifications à apporter seront alors mineures. Cela est plus pertinent que de chercher à utiliser des bibliothèques spécifiquement conçues à cet effet dont la pérennité n'est pas garantie.

In [10]:
# Représentation classique avec un titre par exemple
plot(dep$geometry, col="#EEEEEE", main = "La répartition de la population")
plot(centres, add=TRUE, cex=sqrt(dep$Pop_totale) / 500, pch = 21, bg="#EEA861", lwd = 1, col = "black")

# Créer une légende vide pour réserver l'espace
valeurs_legende <- c(min(dep$Pop_totale), mean(dep$Pop_totale), max(dep$Pop_totale))
leg_info <- legend("bottomleft", 
                   legend = rep("", length(valeurs_legende)), # Texte vide
                   title = "Population",
                   bty = "n",
                   y.intersp = 1.2) # Espace vertical entre les lignes à modifier en fonction du contexte

# Ajouter les cercles sur la légende en utilisant les coordonnées récupérées
tailles_legende <- sqrt(valeurs_legende) / 500
points(rep(leg_info$text$x[1] - 80000 , length(valeurs_legende)), # Position X à modifier en fonction du contexte
       leg_info$text$y,                                       # Position Y
       pch = 21, bg = "#EEA861", col = "black",
       cex = tailles_legende)

# Ajouter les textes de valeurs à côté
text(leg_info$text$x - 40000, leg_info$text$y, 
     labels = round(valeurs_legende, 2), pos = 4)
No description has been provided for this image

On peut aussi s'appuyer sur des points pour réaliser des triangulations, comme celle de Delaunay. On utilisera pour cela la fonction st_triangulate(), son dual les polygones de Voronoi s'obtiennent en utilisant la fonction st_voronoi(). Il est nécessaire pour faire fonctionner ces fonctions de regrouper les points dans un même objet géographique de type multipoint. Pour cela, on utilisera la fonction st_union().

In [11]:
tri <- st_triangulate(st_union(centres))
vor <- st_voronoi(st_union(centres))
In [12]:
plot(dep$geom, col="#EEEEEE")
plot(tri, col=NA, border = "#466EE0",add=TRUE)
plot(centres, add=TRUE)
plot(tri,col=NA,border=NA)
plot(vor, col="#EEA861", add=TRUE)
plot(centres, add=TRUE)
No description has been provided for this image
No description has been provided for this image

Grand classique du traitement de l'information géographique, les buffers peuvent être réalisés en utilisant la fonction st_buffer().

In [13]:
buf <- st_buffer(centres, dist = 20000)
plot(dep$geometry, col="#EEEEEE")
plot(buf, col ="#EEA861", border = "black", add=TRUE)
plot(centres, pch=16, cex = 0.5, add=TRUE, col ="#555555")
No description has been provided for this image

Les intersections s'obtiennent en utilisant la fonction st_intersection(). Vous avez compris, la plupart des fonctions sf ont des noms relativement simples à trouver...

In [14]:
buf2 <- st_buffer(centres, dist = 50000)
inter <- st_intersection(buf2, dep$geometry)
plot(inter)
No description has been provided for this image

Pour les regroupements géométriques, c'est un peu plus subtil. Souvent, en anglais, on utilise le terme dissolve. La bibliothèque sf préfère parler d'union, on utilisera donc la fonction st_union().

In [15]:
france <- st_union(dep$geometry)
plot(france, col ="#EEA861")
No description has been provided for this image

En utilisant des tests pour filtrer ce que l'on souhaite regrouper, on peut réaliser des regroupements plus complexes.

In [16]:
depidf <- dep[dep$CODE_REG=='11',] # Récupération des départements d'Ile-de-France
idf <- st_union(depidf$geometry) # Création d'une région Ile-de-France
plot(idf, col ="#EEA861")
No description has been provided for this image

A l'aide de tests, on peut ainsi réaliser une carte des régions de la France métropolitaine. Néanmoins, cela peut être jugé un peu complexe pour pas grand-chose. On peut alors plutôt utiliser la fonction aggregate(). Cette fonction requiert généralement de se questionner sur comment l'on va gérer l'information sémantique (tabulaire).

In [17]:
regions <- aggregate(
  dep["Pop_totale"],      # On sélectionne la ou les colonnes de données (optionnel)
  by = list(ID = dep$CODE_REG, NOM_REGION = dep$NOM_REGION), #Champs utilisé pour faire le regroupement on peut en utiliser qu'un ici
  FUN = sum                        # Fonction pour traiter les données numériques
)
plot(regions$geometry, col = 2:7)
No description has been provided for this image
In [18]:
print(regions)

Simple feature collection with 22 features and 3 fields
Attribute-geometry relationship: 0 constant, 1 aggregate, 2 identity
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 99227.52 ymin: 6049649 xmax: 1242375 ymax: 7110525
Projected CRS: RGF93 / Lambert-93
First 10 features:
   ID        NOM_REGION Pop_totale                       geometry
1  42            ALSACE    1837087 POLYGON ((1040404 6789657, ...
2  72         AQUITAINE    3177625 MULTIPOLYGON (((427834.5 61...
3  83          AUVERGNE    1341863 POLYGON ((634691.8 6393270,...
4  25   BASSE-NORMANDIE    1467425 MULTIPOLYGON (((533474 6789...
5  26         BOURGOGNE    1638588 POLYGON ((805986.9 6568662,...
6  53          BRETAGNE    3149701 MULTIPOLYGON (((231447.5 67...
7  24            CENTRE    2531588 POLYGON ((577533.3 6585981,...
8  21 CHAMPAGNE-ARDENNE    1338004 POLYGON ((877232 6724190, 8...
9  94             CORSE     302966 MULTIPOLYGON (((1224200 605...
10 43     FRANCHE-COMTE    1163931 MULTIPOLYGON (((934741 6595...

Cette fonction aggregate() est une fonction classique de R. Elle peut s'appliquer à de simples dataFrames, pas besoin d'un format géographique. La bibliothèque sf prend simplement le relais lorsque les objets sont dans un format géographique. Dans le cas classique, donc en dehors des formats géographiques, cette fonction aggregate() a la particularité de fonctionner en utilisant le caractère ~ comme dans la fonction lm().

In [19]:
table_reg <- aggregate(cbind(dep$Pop_hommes, dep$Pop_femmes, dep$Pop_Etrang, dep$Pop_totale) ~ NOM_REGION, data = dep, FUN = sum)
names(table_reg) <- c('NOM_REGION','Pop_hommes','Pop_femmes','Pop_Etrang','Pop_totale')
head(table_reg)
A data.frame: 6 × 5
NOM_REGIONPop_hommesPop_femmesPop_EtrangPop_totale
<chr><dbl><dbl><dbl><dbl>
1ALSACE 897015 9400781394811837087
2AQUITAINE 152801716496291221773177625
3AUVERGNE 649618 692256 436441341863
4BASSE-NORMANDIE 710862 756580 274321467425
5BOURGOGNE 793427 845175 542371638588
6BRETAGNE 15295291620175 628593149701

Avec des données géographiques, si on souhaite minimiser le traitement des données tabulaires pour se concentrer sur les transformations géographiques, la seule information sémantique conservée sera l'identifiant utilisé pour réaliser le regroupement. Ces identifiants pourront servir par la suite pour joindre d'autres données sémantiques disponibles à cette nouvelle échelle. Il faudra alors faire attention à bien choisir cet identifiant.

In [20]:
geom_reg <- aggregate(dep['NOM_REGION'], by = list(NOM_REGION = dep$NOM_REGION), FUN=min)
print(geom_reg)

Simple feature collection with 22 features and 2 fields
Attribute-geometry relationship: 0 constant, 1 aggregate, 1 identity
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 99227.52 ymin: 6049649 xmax: 1242375 ymax: 7110525
Projected CRS: RGF93 / Lambert-93
First 10 features:
          NOM_REGION      NOM_REGION.1                       geometry
1             ALSACE            ALSACE POLYGON ((1040404 6789657, ...
2          AQUITAINE         AQUITAINE MULTIPOLYGON (((427834.5 61...
3           AUVERGNE          AUVERGNE POLYGON ((634691.8 6393270,...
4    BASSE-NORMANDIE   BASSE-NORMANDIE MULTIPOLYGON (((533474 6789...
5          BOURGOGNE         BOURGOGNE POLYGON ((805986.9 6568662,...
6           BRETAGNE          BRETAGNE MULTIPOLYGON (((231447.5 67...
7             CENTRE            CENTRE POLYGON ((577533.3 6585981,...
8  CHAMPAGNE-ARDENNE CHAMPAGNE-ARDENNE POLYGON ((877232 6724190, 8...
9              CORSE             CORSE MULTIPOLYGON (((1224200 605...
10     FRANCHE-COMTE     FRANCHE-COMTE MULTIPOLYGON (((934741 6595...

La fonction st_crs() permet d'obtenir les informations sur la projection de l'objet sf.

In [21]:
 st_crs(dep)

Coordinate Reference System:
  User input: RGF93 / Lambert-93 
  wkt:
PROJCRS["RGF93 / Lambert-93",
    BASEGEOGCRS["RGF93",
        DATUM["Reseau Geodesique Francais 1993",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4171]],
    CONVERSION["Lambert-93",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",46.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",3,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",49,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",44,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",700000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",6600000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["France - onshore and offshore, mainland and Corsica."],
        BBOX[41.15,-9.86,51.56,10.38]],
    ID["EPSG",2154]]

La fonction st_transform() permet de changer de système de coordonnées.

In [22]:
centres_4326  <- st_transform(centres, crs = 4326)
print(centres_4326[1])

Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 5.348801 ymin: 46.09944 xmax: 5.348801 ymax: 46.09944
Geodetic CRS:  WGS 84

POINT (5.348801 46.09944)

Sélection et jointure¶

On a fait un bon tour des principaux traitements géométriques. Désormais, on peut s'intéresser à l'enrichissement des données en s'intéressant par exemple aux jointures attributaires. La fonction classiquement utilisée pour cela est la fonction merge(). L'argument by permet de définir la colonne commune permettant d'effectuer la jointure.

In [23]:
resultat <- merge(table_reg, geom_reg, by = "NOM_REGION")
head(resultat)
A data.frame: 6 × 7
NOM_REGIONPop_hommesPop_femmesPop_EtrangPop_totaleNOM_REGION.1geometry
<chr><dbl><dbl><dbl><dbl><chr><GEOMETRY [m]>
1ALSACE 897015 9400781394811837087ALSACE POLYGON ((1040404 6789657, ...
2AQUITAINE 152801716496291221773177625AQUITAINE MULTIPOLYGON (((427834.5 61...
3AUVERGNE 649618 692256 436441341863AUVERGNE POLYGON ((634691.8 6393270,...
4BASSE-NORMANDIE 710862 756580 274321467425BASSE-NORMANDIEMULTIPOLYGON (((533474 6789...
5BOURGOGNE 793427 845175 542371638588BOURGOGNE POLYGON ((805986.9 6568662,...
6BRETAGNE 15295291620175 628593149701BRETAGNE MULTIPOLYGON (((231447.5 67...

Si les données ne proviennent pas de la même source, il est peu probable que les noms des champs contenant les identifiants permettant d'effectuer la jointure attributaire soient les mêmes. On peut ainsi charger les données des catégories socio-professionnelles disponibles ici : https://sergelhomme.fr/data/TPR.csv

In [24]:
ad <- file.choose()
csp <- read.csv(ad, sep=';')
head(csp)
A data.frame: 6 × 10
IDNOMAgriculteurArtisan.CommercantCadreProfession.IntermediaireEmployeOuvrierRetraiteAutres
<chr><chr><int><int><int><int><int><int><int><int>
101Ain 40671674536426 70663 77349789983041729910
202Aisne 52011141418629 48889 71578789063206342108
303Allier 6159 939611842 31102 46196421012421821344
404Alpes-de-Haute-Provence2027 6387 6951 16140 20885155601068710278
505Hautes-Alpes 2040 5286 5883 15687 1970712415 8865 7222
606Alpes-Maritimes 19183951473884117652160574872275480168990

Les arguments by.x et by.y de la fonction merge() permettent alors de définir les noms des champs à utiliser pour les deux sources de données afin de réaliser la jointure.

In [25]:
dep_csp <- merge(dep, csp, by.x = "CODE_DEPT", by.y = "ID")
print(dep_csp)

Simple feature collection with 96 features and 26 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 99227.52 ymin: 6049649 xmax: 1242375 ymax: 7110525
Projected CRS: RGF93 / Lambert-93
First 10 features:
   CODE_DEPT ID_GEOFLA                NOM_DEPT              NOM_CHF CODE_REG
1         01         1                     AIN      BOURG-EN-BRESSE       82
2         02         2                   AISNE                 LAON       22
3         03         3                  ALLIER              MOULINS       83
4         04         4 ALPES-DE-HAUTE-PROVENCE      DIGNE-LES-BAINS       93
5         05         5            HAUTES-ALPES                  GAP       93
6         06         6         ALPES-MARITIMES                 NICE       93
7         07         7                 ARDECHE               PRIVAS       82
8         08         8                ARDENNES CHARLEVILLE-MEZIERES       21
9         09         9                  ARIEGE                 FOIX       73
10        10        10                    AUBE               TROYES       21
                   NOM_REGION Pop_hommes Pop_femmes Pop_Etrang Pop_totale
1                 RHONE-ALPES     287336     294023      45109     581355
2                    PICARDIE     262006     276794      13447     538790
3                    AUVERGNE     164392     178416       9394     342807
4  PROVENCE-ALPES-COTE-D'AZUR      76580      81387       6388     157965
5  PROVENCE-ALPES-COTE-D'AZUR      65494      68714       4192     134205
6  PROVENCE-ALPES-COTE-D'AZUR     511716     572712     100388    1084428
7                 RHONE-ALPES     151678     159776       9343     311452
8           CHAMPAGNE-ARDENNE     139346     144856      11382     284197
9               MIDI-PYRENEES      73313      76894       6864     150201
10          CHAMPAGNE-ARDENNE     146609     154726      14327     301327
   Naissances  Deces Solde_nat Solde_mig Tx_Etrange Superfice      Sup
1       61000  37114     23886     41991     0.0776  5784.556 5773.960
2       62669  48087     14582    -11105     0.0250  7419.598 7426.966
3       30652  38838     -8186      6378     0.0274  7380.145 7366.255
4       13741  13972      -231     18513     0.0404  6995.771 6994.779
5       12882  10793      2089     10485     0.0312  5689.465 5684.619
6      102510 101409      1101     71461     0.0926  4295.464 4295.917
7       29910  28091      1819     23473     0.0300  5565.608 5560.374
8       32094  25080      7014    -12941     0.0400  5243.878 5249.568
9       12445  15900     -3455     16309     0.0457  4906.472 4916.098
10      31923  25772      6151      3077     0.0475  6025.735 6020.256
                       NOM Agriculteur Artisan.Commercant Cadre
1                      Ain        4067              16745 36426
2                    Aisne        5201              11414 18629
3                   Allier        6159               9396 11842
4  Alpes-de-Haute-Provence        2027               6387  6951
5             Hautes-Alpes        2040               5286  5883
6          Alpes-Maritimes        1918              39514 73884
7                  Ardeche        4132              10758 12461
8                 Ardennes        3362               6520  9771
9                   Ariege        2348               4888  5537
10                    Aube        4786               7198 12499
   Profession.Intermediaire Employe Ouvrier Retraite Autres
1                     70663   77349   78998    30417  29910
2                     48889   71578   78906    32063  42108
3                     31102   46196   42101    24218  21344
4                     16140   20885   15560    10687  10278
5                     15687   19707   12415     8865   7222
6                    117652  160574   87227    54801  68990
7                     31859   39088   37199    21660  19063
8                     25487   36326   42618    16464  23226
9                     14442   20788   16635    10379   9741
10                    28889   40121   43224    18616  18373
                         geometry
1  MULTIPOLYGON (((919195 6541...
2  MULTIPOLYGON (((735604.5 68...
3  MULTIPOLYGON (((753768.9 65...
4  MULTIPOLYGON (((992637.8 63...
5  MULTIPOLYGON (((1012912 640...
6  MULTIPOLYGON (((1018256 627...
7  MULTIPOLYGON (((831641.3 63...
8  MULTIPOLYGON (((842093.3 69...
9  MULTIPOLYGON (((631546.3 61...
10 MULTIPOLYGON (((796594.5 67...

Avec cette jointure, on peut alors cartographier les taux d'ouvriers par département. Pour réaliser une cartographie en aplat de couleurs, on peut utiliser la bibliothèque RColorBrewer et sa fonction brewer.pal() pour obtenir les codes hexadécimaux des couleurs et la bibliothèque classInt et sa fonction classIntervals() pour obtenir les seuils des classes de la discrétisation souhaitée. Pour rappel, pour obtenir une liste de couleurs pouvant être utilisée dans un plot R à partir des choix de palette et de discrétisation, il faut utiliser la cut() combinée à la fonction as.character() (revoir le notebook précédent).

In [26]:
Tx_Ouv <- dep_csp$Ouvrier / (dep_csp$Ouvrier + dep_csp$Cadre + dep_csp$Employe + dep_csp$Agriculteur + dep_csp$Artisan.Commercant + dep_csp$Profession.Intermediaire + dep_csp$Autres + dep_csp$Retraite)
library(RColorBrewer)
library(classInt)
palette <- brewer.pal(n = 5, name = "Purples")
br <- classIntervals(Tx_Ouv, n=5, style="jenks")
a <- as.character(cut(Tx_Ouv, breaks = br$brks, labels = palette))
plot(dep_csp$geometry, col = a, border = "black", main="Les Taux d'ouvrier.e.s")
No description has been provided for this image

Pour gérer la légende de cet aplat de couleurs, ce sera plus simple que pour les symboles proportionnels, car on dispose d'une fonction findColours() dont les résultats peuvent être exploités dans la fonction legend().

In [27]:
br$brks <- round(br$brks, 2) # Arrondi des seuils pour éviter d'avoir une légende illisible
col_points <- findColours(br, palette)
plot(dep_csp$geometry, col = a, border = "black", main="Répartition des ouvrier.e.s")
legend(
  "bottomleft", 
  legend = names(attr(col_points, "table")), # Récupère les labels des classes
  fill = attr(col_points, "palette"),       # Récupère les couleurs de la palette
  title = "Taux Ouvrier.e",
  bty = "n"                                 # Supprime le cadre de la légende
)
No description has been provided for this image

Pour sélectionner (identifier) un ou des objets géographiques à l'aide de leurs attributs, on peut utiliser les tests comme on le ferait avec un dataFrame classique (revoir le notebook précédent). Pour les objets géographiques de la classe sf, ces tests peuvent être utilisés dans la fonction plot() de telle manière à identifier cartographiquement des objets géographiques.

In [28]:
print(dep[dep$CODE_DEPT == '75',])

Simple feature collection with 1 feature and 17 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 643057.3 ymin: 6857473 xmax: 661070.4 ymax: 6867035
Projected CRS: RGF93 / Lambert-93
   ID_GEOFLA CODE_DEPT NOM_DEPT                   NOM_CHF CODE_REG
76        76        75    PARIS PARIS--1ER-ARRONDISSEMENT       11
      NOM_REGION Pop_hommes Pop_femmes Pop_Etrang Pop_totale Naissances  Deces
76 ILE-DE-FRANCE    1042203    1169094     329922    2211297     284597 138230
   Solde_nat Solde_mig Tx_Etrange Superfice                       geometry
76    146367    -60921     0.1492   105.365 MULTIPOLYGON (((657172.4 68...
        Sup
76 105.3418
In [29]:
print(dep[dep$Pop_totale > 1000000,c('NOM_DEPT','Pop_totale')])

Simple feature collection with 22 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 281034.3 ymin: 6179669 xmax: 1082669 ymax: 7110525
Projected CRS: RGF93 / Lambert-93
First 10 features:
           NOM_DEPT Pop_totale                       geometry
6   ALPES-MARITIMES    1084428 MULTIPOLYGON (((1018256 627...
13 BOUCHES-DU-RHONE    1966005 MULTIPOLYGON (((917346.8 62...
32    HAUTE-GARONNE    1217344 MULTIPOLYGON (((524704.2 61...
34          GIRONDE    1421276 MULTIPOLYGON (((449162.6 63...
35          HERAULT    1019798 MULTIPOLYGON (((719561.6 62...
39            ISERE    1188660 MULTIPOLYGON (((921862.8 64...
45 LOIRE-ATLANTIQUE    1255871 MULTIPOLYGON (((385081 6667...
58          MOSELLE    1042230 MULTIPOLYGON (((1038623 689...
60             NORD    2564959 MULTIPOLYGON (((781907.7 69...
63    PAS-DE-CALAIS    1459531 MULTIPOLYGON (((706472.2 69...
In [30]:
plot(dep[dep$Pop_totale <= 1000000,]$geom, col="#DDDDDD")
plot(dep[dep$Pop_totale > 1000000,]$geom, col="#EEA861", add=TRUE)
No description has been provided for this image

Pour les sélections spatiales, il y a autant de fonctions que de critères spatiaux. On a ainsi des fonctions st_intersects(), st_contains(), st_covered_by(), st_crosses(), st_equals(), st_is_within_distance(), st_nearest_feature(), st_overlaps(), st_touches(), st_within()...

In [31]:
index <- st_intersects(centres, depidf)
print(index)
print(index[76]) # Pour Paris

Sparse geometry binary predicate list of length 96, where the predicate
was `intersects'
first 10 elements:
 1: (empty)
 2: (empty)
 3: (empty)
 4: (empty)
 5: (empty)
 6: (empty)
 7: (empty)
 8: (empty)
 9: (empty)
 10: (empty)
[[1]]
[1] 1

Comme vous le voyez, par défaut, ces fonctions renvoient une liste qui contient soit des éléments vides pour les éléments de la couche 1 qui ne sont pas concernés par la relation spatiale avec la couche 2, soit les index des éléments de la couche 2 qui correspondent à la relation spatiale choisie. Il peut être plus simplement utile de savoir si oui ou non les éléments de la couche 1 sont concernés par la relation spatiale étudiée avec les éléments de la couche 2. Pour cela, on peut utiliser l'argument sparse.

In [32]:
index <- st_intersects(centres, depidf, sparse=FALSE)
plot(idf, col="#EEA861")
plot(centres[rowSums(index) > 0], pch=21, add=TRUE, bg='grey')
No description has been provided for this image

Pour effectuer des jointures spatiales, on utilisera la fonction st_join(). L'argument join de cette fonction permet de définir la relation spatiale souhaitée. Ces relations spatiales correspondent au nom des fonctions permettant d'effectuer les sélections spatiales correspondantes.

In [33]:
pt_csp <- st_join(st_sf(centres), dep_csp, join = st_intersects)
print(pt_csp[,c('NOM_DEPT','Employe')])

Simple feature collection with 96 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 176848 ymin: 6103614 xmax: 1211296 ymax: 7044252
Projected CRS: RGF93 / Lambert-93
First 10 features:
                  NOM_DEPT Employe                 geometry
1                      AIN   77349 POINT (881425.4 6558215)
2                    AISNE   71578   POINT (740412 6940173)
3                   ALLIER   46196 POINT (714474.9 6588217)
4  ALPES-DE-HAUTE-PROVENCE   20885 POINT (959609.2 6339447)
5             HAUTES-ALPES   19707 POINT (958583.4 6401404)
6          ALPES-MARITIMES  160574  POINT (1030360 6324001)
7                  ARDECHE   39088 POINT (812764.2 6406858)
8                 ARDENNES   36326 POINT (818622.2 6947516)
9                   ARIEGE   20788 POINT (577729.8 6203552)
10                    AUBE   40121 POINT (786149.2 6801093)

Spdep¶

Spdep est la bibliothèque R classiquement utilisée pour étudier les dépendances spatiales entre objets géographiques. Elle permet notamment de mesurer de l'autocorrélation spatiale.

In [ ]:
library(spdep)

Pour mesurer de l'autocorrélation spatiale, il faut définir un voisinage spatial. Pour cela, on peut utiliser la fonction poly2nb() qui permet de définir des relations de contiguités pour des polygones.

In [35]:
depqueen <- poly2nb(dep$geometry)

L'affichage de ces relations spatiales passent par la création de points à partir desquels tracer les relations. On peut s'appuyer pour cela sur les centroides.

In [36]:
coords <- st_centroid(dep$geometry)
plot(dep$geometry)
plot(depqueen,coords, col = 'red',add=TRUE)
No description has been provided for this image

Ensuite, la fonction nb2listw() transforme des relations de voisinage en une pondération spatiale que l'on peut utiliser pour le calcul d'autocorrélation spatiale qui s'effectuera avec la fonction moran.test(). Ici, on utilise les taux d'ouvriers calculés plus haut dans ce notebook pour effectuer le calcul.

In [37]:
depqueenw <- nb2listw(depqueen)
Iqueen <- moran.test(Tx_Ouv,depqueenw)
print(Iqueen)

	Moran I test under randomisation

data:  Tx_Ouv  
weights: depqueenw    

Moran I statistic standard deviate = 7.9693, p-value = 7.977e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.52632636       -0.01052632        0.00453802

La fonction moran.plot() permet d'afficher le diagramme de Moran correspondant.

In [38]:
moran.plot(Tx_Ouv,depqueenw)
No description has been provided for this image

Spdep est une bibliothèque très connue, car elle permet de définir de nombreux types de voisinage. En effet, la bibliothèque Spdep n'est pas uniquement utilisée pour calculer de l'autocorrélation spatiale. Ses fonctions de définition de voisinage et de pondération spatiale peuvent être utilisées pour d'autres calculs.

In [39]:
iristri <- tri2nb(coords) # Triangulation de Delaunay
plot(iristri, coords, col = 'red')
title(main = "Voisinage par triangulation de Delaunay")
irisdn <- dnearneigh(coords,0, 100000) # Distance maximale
plot(irisdn,coords, col = 'red')
title(main = "Voisinage par distance maximale (100km)")
knn <- knearneigh(coords, k=1) # Plus proche voisin
irisknn <- knn2nb(knn)
plot(irisknn,coords,col='red')
title(main = "Voisinage au plus proche voisin")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

La fonction localmoran() permet de calculer des indices de Moran locaux.

In [40]:
lisamoran <- localmoran(Tx_Ouv, depqueenw)
In [41]:
palette <- brewer.pal(n=5,name = "RdBu")
palette <- rev(palette)
br <- c(min(lisamoran[,1]), -0.4, -0.2, 0.2, 0.4, max(lisamoran[,1]) )
couleur <- as.character(cut(lisamoran[,1], labels=palette, breaks=br))
plot(dep$geometry,col=couleur)
No description has been provided for this image

On peut alors identifier des hotspots et des coldspots.

In [42]:
tx_ouv_voi <- lag(depqueenw, Tx_Ouv)
Tx_Ouv_std <- scale(Tx_Ouv)
tx_ouv_voi_std <- scale(tx_ouv_voi)
In [43]:
couleur <- rep("#EEEEEE",length(tx_ouv_voi))
couleur[lisamoran[,5]<0.05& Tx_Ouv_std < 0 & tx_ouv_voi_std < 0] <- "#0000FF"
couleur[lisamoran[,5]<0.05 & Tx_Ouv_std > 0 & tx_ouv_voi_std > 0] <- "#FF0000"
plot(dep$geometry, col=couleur)
No description has been provided for this image
In [44]:
plot(dep$geometry, col=couleur, main = "Les Ouvrier.e.s en France métropolitaine")
legend(
  "bottomleft", 
  legend = c("HotSpot","ColdSpot"), # Récupère les labels des classes
  fill = c("#FF0000","#0000FF"),       # Récupère les couleurs de la palette
  title = "Taux Ouvrier.e",
  bty = "n"                                 # Supprime le cadre de la légende
)
No description has been provided for this image

Sp¶

Enfin, la bibliothèque Sp est une bibliothèque classique pour traiter de l'information géographique, car elle propose un type (une classe) d'objets R très utilisé par des fonctions R.

In [45]:
library(sp)

La fonction as_spatial() permet alors de transformer un objet sf en un objet sp. L'affichage d'un objet sp se fait alors plus simplement dans un plot, pas besoin d'aller chercher une colonne géométrie comme pour sf.

In [46]:
carte <- as_Spatial(dep)
plot(carte)
No description has been provided for this image

La classe d'objets sp repose sur des slots que l'on peut aller chercher avec le caractère @ . Ainsi pour avoir accès aux tables attributaires, on utilise @data, puis on va chercher les champs avec les $ comme pour un dataFrame.

In [47]:
carte@data$Pop_totale
  1. 581355
  2. 538790
  3. 342807
  4. 157965
  5. 134205
  6. 1084428
  7. 311452
  8. 284197
  9. 150201
  10. 301327
  11. 349237
  12. 275889
  13. 1966005
  14. 678206
  15. 148737
  16. 351581
  17. 611714
  18. 313251
  19. 242896
  20. 521608
  21. 581570
  22. 123907
  23. 409388
  24. 522685
  25. 478069
  26. 577087
  27. 423559
  28. 890509
  29. 140953
  30. 162013
  31. 694323
  32. 1217344
  33. 185266
  34. 1421276
  35. 1019798
  36. 967588
  37. 232004
  38. 585406
  39. 1188660
  40. 260740
  41. 373142
  42. 326599
  43. 742076
  44. 221834
  45. 1255871
  46. 650769
  47. 172796
  48. 326399
  49. 76973
  50. 774823
  51. 496937
  52. 566010
  53. 186470
  54. 302983
  55. 729768
  56. 194218
  57. 710034
  58. 1042230
  59. 220653
  60. 2564959
  61. 799725
  62. 292282
  63. 1459531
  64. 628485
  65. 647420
  66. 229079
  67. 441387
  68. 1091015
  69. 746072
  70. 1690498
  71. 238548
  72. 553968
  73. 559587
  74. 408842
  75. 716277
  76. 2211297
  77. 1248580
  78. 1303702
  79. 1406053
  80. 365059
  81. 568086
  82. 371738
  83. 235915
  84. 1001408
  85. 538902
  86. 616906
  87. 424354
  88. 373940
  89. 380145
  90. 342359
  91. 141958
  92. 1205850
  93. 1549619
  94. 1506466
  95. 1310876
  96. 1165397

Néanmoins, désormais, on peut directement utiliser le $ pour avoir accès aux champs de données.

In [48]:
carte$Pop_totale
  1. 581355
  2. 538790
  3. 342807
  4. 157965
  5. 134205
  6. 1084428
  7. 311452
  8. 284197
  9. 150201
  10. 301327
  11. 349237
  12. 275889
  13. 1966005
  14. 678206
  15. 148737
  16. 351581
  17. 611714
  18. 313251
  19. 242896
  20. 521608
  21. 581570
  22. 123907
  23. 409388
  24. 522685
  25. 478069
  26. 577087
  27. 423559
  28. 890509
  29. 140953
  30. 162013
  31. 694323
  32. 1217344
  33. 185266
  34. 1421276
  35. 1019798
  36. 967588
  37. 232004
  38. 585406
  39. 1188660
  40. 260740
  41. 373142
  42. 326599
  43. 742076
  44. 221834
  45. 1255871
  46. 650769
  47. 172796
  48. 326399
  49. 76973
  50. 774823
  51. 496937
  52. 566010
  53. 186470
  54. 302983
  55. 729768
  56. 194218
  57. 710034
  58. 1042230
  59. 220653
  60. 2564959
  61. 799725
  62. 292282
  63. 1459531
  64. 628485
  65. 647420
  66. 229079
  67. 441387
  68. 1091015
  69. 746072
  70. 1690498
  71. 238548
  72. 553968
  73. 559587
  74. 408842
  75. 716277
  76. 2211297
  77. 1248580
  78. 1303702
  79. 1406053
  80. 365059
  81. 568086
  82. 371738
  83. 235915
  84. 1001408
  85. 538902
  86. 616906
  87. 424354
  88. 373940
  89. 380145
  90. 342359
  91. 141958
  92. 1205850
  93. 1549619
  94. 1506466
  95. 1310876
  96. 1165397

La classe sp propose des mesures d'objets très appréciées, comme le calcul des superficies pour les polygones. En manipulant les slots de cette classe, on peut alors se passer d'effectuer certains calculs supplémentaires.

In [49]:
carte@polygons[[1]]@area / 1000000
5773.95987277597