voisinage de polygone et attributs

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Mathieu Lagarde
Messages : 10
Enregistré le : 19 Mar 2010, 09:28

voisinage de polygone et attributs

Messagepar Mathieu Lagarde » 07 Sep 2020, 12:22

Bonjour,
Je sèche sur un problème de voisinage de polygones...
Explications...
Je pars d'un shapefile avec les départements français (polygones) auxquels j'ai ajouté des variables / colonnes. Chacune de ces variables est une espèce codée en "0" (absence dans le département) ou "1" (présence dans le département).

Code : Tout sélectionner

> rt <- readOGR("hist_rt.shp",layer="hist_rt")
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\mlagarde\Documents\R\hist_rt.shp", layer: "hist_rt"
with 96 features
It has 151 fields
Integer64 fields read as strings:  ### liste des taxons ###


Code : Tout sélectionner

> class(rt)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"


Code : Tout sélectionner

> head(rt)
class       : SpatialPolygonsDataFrame
features    : 6
extent      : 644461.3, 1077719, 6272508, 6997025  (xmin, xmax, ymin, ymax)
crs         : +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs
variables   : 151
names       :                       ID,      NOM_DEP,    NOM_DEP_M, INSEE_DEP, INSEE_REG, CHF_DEP,      DEP_NOM, tax1, tax2, tax3, tax4, tax5, tax6, tax7, tax8, ...
min values  : DEP000000000000000000001,          Ain,          AIN,        01,        32,   01053,          Ain,        0,        0,        1,        0,        0,        0,        0,        1, ...
max values  : DEP000000000000000000006, Hautes-Alpes, HAUTES ALPES,        06,        93,   06088, Hautes-Alpes,        2,        2,        2,        2,        2,        2,        2,        2, ...



Mon objectif est d'ajouter un paramètre à ces variables: "2" (présence potentielle dans le département), selon une règle très simple (et discutable): si le département A est codé "0" et est limitrophe à un département codé "1", alors le département A doit être codé "2".

Je me suis penché sur la fonction poly2nb du package spdep qui me semble répondre à mes besoins, jusqu'à obtenir la matrice de voisinage des départements:

Code : Tout sélectionner

> head(mat)
                        Ain Aisne Allier Alpes-de-Haute-Provence Hautes-Alpes Alpes-Maritimes Ardèche Ardennes Ariège
Ain                       0     0      0                       0            0               0        0        0       0
Aisne                     0     0      0                       0            0               0        0        1       0
Allier                    0     0      0                       0            0               0        0        0       0
Alpes-de-Haute-Provence   0     0      0                       0            1               1        0        0       0
Hautes-Alpes              0     0      0                       1            0               0        0        0       0
Alpes-Maritimes           0     0      0                       1            0               0        0        0       0


Et c'est là que je sèche. Comment puis-je exploiter cette matrice pour attribuer le nouveau paramètre à chacune de mes variables "espèces"?
Et avant tout, est-ce que j'ai fait le bon choix de partir sur poly2nb ?

D'avance, merci pour vos pistes !

Sébastien Rochette
Messages : 8
Enregistré le : 03 Juil 2020, 12:43

Re: voisinage de polygone et attributs

Messagepar Sébastien Rochette » 07 Sep 2020, 12:56

Bonjour,

La matrice donne le résultat recherché, donc a priori, c'est une fonction adaptée.
En revanche, pour la lecture d'un shapefile, je vous recommande d'utiliser le package {sf}. Non seulement, {sp} est à considérer comme obsolète, mais en plus, les données seront plus faciles à manipuler. Vous travailler avec les données comme sur une table de données classique directement.

Pour votre exemple, je ferai la somme des lignes de la matrice et j'en ferai une colonne d'une nouvelle table de données. Cette nouvelle table contient une colonne "nom_departement" et une colonne "somme_des_lignes".
Ensuite, vous faites une jointure entre les données carto et cette table.
A partir de là, c'est de la comparaison de valeur dans une table de données classique. (A condition d'utiliser {sf})

Pour que je puisse vous apporter plus d'aide, il me faudrait une exemple reproductible, avec les données disponibles. cf: https://thinkr.fr/reprex-ou-comment-demander-de-laide-efficacement/

Serge Rapenne
Messages : 1386
Enregistré le : 20 Aoû 2007, 15:17
Contact :

Re: voisinage de polygone et attributs

Messagepar Serge Rapenne » 07 Sep 2020, 14:34

Bonjour,

je pense que la fonction "gtouches" du package "rgeos" devrait te permettre de faire ça. Sans les données, il est effectivement difficile de te donner du code mais, la logique est de parcourir les départements un à un, de trouver les limitrophes avec "gtouches", de changer la valeur des variables si besoin puis de passer au suivant.

Serge

Mathieu Lagarde
Messages : 10
Enregistré le : 19 Mar 2010, 09:28

Re: voisinage de polygone et attributs

Messagepar Mathieu Lagarde » 10 Sep 2020, 09:46

Bonjour,
Merci pour vos réponses.
Afin d'obtenir un reprex, j'ai revu mes données de base et suis parti du package {raster} pour obtenir les départements français.
Voici le code qui me permet d'arriver au même point que mon message de départ.

Code : Tout sélectionner

# chargement des packages
library(raster)
library(spdep)
# lecture des données
fr <- getData(name = "GADM", country = "FRA", level = 2) # carte des département français
fridf <- subset(fr, NAME_1 == "Île-de-France")

data_fict <- tibble::tribble(
  ~CC_2, ~tax1, ~tax2, ~tax3,
   "75",     1,     0,     0,
   "77",     1,     0,     1,
   "78",     0,     0,     1,
   "91",     0,     0,     0,
   "92",     0,     1,     0,
   "93",     0,     1,     0,
   "94",     0,     1,     0,
   "95",     0,     0,     0
)

# jointure attributaire
jfrdat <- merge(fridf, data_fict, by = "CC_2")

# utilisation de la fonction poly2nb pour obtenir la liste des voisins
rt_nb <- poly2nb(jfrdat,
  row.names = jfrdat$CC_2, snap = sqrt(.Machine$double.eps),
  queen = TRUE, useC = TRUE, foundInBox = NULL
)
mat <- nb2mat(rt_nb, style = "B")
colnames(mat) <- rownames(mat)
mat


Par contre, je ne comprends pas où ce que vous me proposez me permettra d'arriver à mes fins:

Sébastien Rochette a écrit :Pour votre exemple, je ferai la somme des lignes de la matrice et j'en ferai une colonne d'une nouvelle table de données. Cette nouvelle table contient une colonne "nom_departement" et une colonne "somme_des_lignes".
Ensuite, vous faites une jointure entre les données carto et cette table.
A partir de là, c'est de la comparaison de valeur dans une table de données classique.


Ce que je souhaite obtenir in fine, c'est cette table:

Code : Tout sélectionner

data_attendue <- tibble::tribble(
  ~CC_2, ~tax1, ~tax2, ~tax3,
   "75",     1,     2,     0,
   "77",     1,     2,     1,
   "78",     0,     2,     1,
   "91",     2,     2,     2,
   "92",     2,     1,     2,
   "93",     2,     1,     2,
   "94",     2,     1,     2,
   "95",     0,     2,     2
)

Sébastien Rochette
Messages : 8
Enregistré le : 03 Juil 2020, 12:43

Re: voisinage de polygone et attributs

Messagepar Sébastien Rochette » 11 Sep 2020, 08:33

Merci pour le reprex, c'est un peu plus facile de travailler avec.
Je vous propose une solution. Toutes les étapes sont décrites. Je vous recommande de tester en ajoutant chaque étape une par une pour voir ce qu'il se passe.

Code : Tout sélectionner

# chargement des packages
library(raster)
library(spdep)
# lecture des données
fr <- getData(name = "GADM", country = "FRA", level = 2) # carte des département français
fridf <- subset(fr, NAME_1 == "Île-de-France")

data_fict <- tibble::tribble(
  ~CC_2, ~tax1, ~tax2, ~tax3,
  "75",     1,     0,     0,
  "77",     1,     0,     1,
  "78",     0,     0,     1,
  "91",     0,     0,     0,
  "92",     0,     1,     0,
  "93",     0,     1,     0,
  "94",     0,     1,     0,
  "95",     0,     0,     0
)

# jointure attributaire
jfrdat <- merge(fridf, data_fict, by = "CC_2")

# utilisation de la fonction poly2nb pour obtenir la liste des voisins
rt_nb <- poly2nb(jfrdat,
                 row.names = jfrdat$CC_2, snap = sqrt(.Machine$double.eps),
                 queen = TRUE, useC = TRUE, foundInBox = NULL
)
mat <- nb2mat(rt_nb, style = "B")
colnames(mat) <- rownames(mat)
mat

# Prenons mat[1,], il dit que 91 est limitrophe à 92, 77, 94, 78
# Je tranforme la matrice en table de données tidy

library(dplyr)
library(tidyr)
library(sf)

# Transform jfrdat en objet {sf} plutôt que {sp} pour lisibilité et utilisation avec {dplyr}
jfrdat_sf <- st_as_sf(jfrdat)
# Transformation de l'objet {sf} en tibble et réduction des variables pour accélérer les opérations
# On rapprochera plus tard de la géométrie
jfrdat_tbl <- st_drop_geometry(jfrdat_sf) %>%
  # utiliser dplyr::select sinon conflit avec raster::select
  dplyr::select("CC_2", starts_with("tax"))

# On combine les variables
tax_combine <-
# Matrice au format tibble
bind_cols(
  # _Dpt focus
  tibble(origin = colnames(mat)),
  # _voisins
  as_tibble(mat)
) %>%
  # Transformation au format long
  pivot_longer(cols = -origin, names_to = "dpt", values_to = "is_neighbor") %>%
  # Filtrer seulement les dpt qui sont des voisins, les autres ne sont pas utiles
  filter(is_neighbor == 1) %>%
  # Ajout des valeurs originales de tax1, tax2, tax3 des dpt voisins par jointure
  # avec "jfrdat_tbl" sur la variable "dpt" (celui du voisin)
  left_join(
    jfrdat_tbl, by = c("dpt" = "CC_2")
  ) %>%
  # Pour chaque origin et chaque tax, on teste s'il y a au moins un 1 chez les voisins
  group_by(origin) %>%
  summarise(
    tax1_ngb = sum(tax1) >= 1,
    tax2_ngb = sum(tax2) >= 1,
    tax3_ngb = sum(tax3) >= 1
  ) %>%
  # J'ajoute les valeurs de "tax" pour le departement d'origine, pour comparer à ceux des voisins
  # Par jointure avec "jfrdat_tbl" mais sur la variable "origin" ctte fois
  left_join(
    jfrdat_tbl, by = c("origin" = "CC_2")
  ) %>%
  # On peut comparer les valeurs de tax une à une
  # e.g. Si tax1 == 1, alors tax1 reste 1
  # Si tax1 == 0 et tax1_ngb == FALSE, alors tax1 reste 0
  # Si tax1 == 0 et tax1_ngb == TRUE, alors devient 2
  # Pour la clarté de l'opération, on crée une nouvelle variable
  mutate(
    tax1_combine = if_else(tax1 == 0 & tax1_ngb == TRUE, 2, tax1),
    tax2_combine = if_else(tax2 == 0 & tax2_ngb == TRUE, 2, tax2),
    tax3_combine = if_else(tax3 == 0 & tax3_ngb == TRUE, 2, tax3)
  )


# On peut maintenant réinjecter le tout dans le jeu de données spatial par jointure

jfrdat_sf_combine <- jfrdat_sf %>%
  left_join(tax_combine, by = c("CC_2" = "origin")) %>%
  # On transforme en character pour la carte
  mutate_at(vars(starts_with("tax")), as.character)

# Et faire les graphs
library(ggplot2)

ggplot(jfrdat_sf_combine) +
  geom_sf(aes(fill = tax1_combine)) +
  scale_fill_viridis_d()


Image

Facundo Muñoz
Messages : 65
Enregistré le : 04 Juil 2019, 09:58
Contact :

Re: voisinage de polygone et attributs

Messagepar Facundo Muñoz » 11 Sep 2020, 08:47

Bonjour,

et bien, le temps que j'ai pris à trouver une solution, Sébastien a répondu.
Tant pis, pas grave. Je vous laisse ma solution quand même car elle est différente. J'utilise seulement les packages de base et ceux qui vous avez utilisé. Vous prendrez ce qui vous arrange.

J'ai juste programmé deux petites fonctions : presence_in_any_neighbour, qui donne un vecteur logique qu'indique la présence d'une espèce dans le voicinage, et recode_species, qui applique votre règle de re-codage à toutes les variables séléctionnées.

Code : Tout sélectionner

# chargement des packages
library(raster)
library(spdep)
# lecture des données
fr <- getData(name = "GADM", country = "FRA", level = 2) # carte des département français
fridf <- subset(fr, NAME_1 == "Île-de-France")

data_fict <- tibble::tribble(
  ~CC_2, ~tax1, ~tax2, ~tax3,
  "75",     1,     0,     0,
  "77",     1,     0,     1,
  "78",     0,     0,     1,
  "91",     0,     0,     0,
  "92",     0,     1,     0,
  "93",     0,     1,     0,
  "94",     0,     1,     0,
  "95",     0,     0,     0
)

#' jointure attributaire
jfrdat <- merge(fridf, data_fict, by = "CC_2")

#' utilisation de la fonction poly2nb pour obtenir la liste des voisins
rt_nb <- poly2nb(jfrdat,
                 row.names = jfrdat$CC_2, snap = sqrt(.Machine$double.eps),
                 queen = TRUE, useC = TRUE, foundInBox = NULL
)

#' Présence d’une espèce dans le voicinage
#'
#' @param x Integer vector. 0 ou 1: absence ou présence dans chaque unité
#' @param n Objet de classe “nb”. Structure de voisinage.
#'
#' @return Veteur logique qui indique la présence de l’espèce dans le voisinage
#' de chaque unité.
presence_in_any_neighbour <- function(x, nb) {
  vapply(nb, function(.) any(x[.] == 1L), TRUE)
}

#' Recodage d’un SpatialPolygonsDataFrame
recode_species <- function(x, vars, nb) {
 
  ans <- x
  for (v in vars) {
    idx <- x[[v]] == 0L & presence_in_any_neighbour(x[[v]], nb)
    ans[[v]][idx] <- 2
  }
 
  return(ans)
}

taxa <- c("tax1", "tax2", "tax3")
jfr_recoded <- recode_species(jfrdat, vars = taxa, nb = rt_nb)

# Vérification du résultat
cbind(
  slot(jfrdat[,c("CC_2", taxa)], "data"),
  slot(jfr_recoded[, taxa], "data")
)[order(jfrdat$CC_2),]
#>   CC_2 tax1 tax2 tax3 tax1 tax2 tax3
#> 1   75    1    0    0    1    2    0
#> 2   77    1    0    1    1    2    1
#> 3   78    0    0    1    0    2    1
#> 4   91    0    0    0    2    2    2
#> 5   92    0    1    0    2    1    2
#> 6   93    0    1    0    2    1    2
#> 7   94    0    1    0    2    1    2
#> 8   95    0    0    0    2    2    2


Vous verrez que ça donne ce que vous attendiez, sauf pour l'unité 95 pour la tax1. Mais je crois que c'est un erreur de votre part, car 95 est voisine de 77 qui continent tax1.

Cordialement,
ƒacu.-

Mathieu Lagarde
Messages : 10
Enregistré le : 19 Mar 2010, 09:28

Re: voisinage de polygone et attributs

Messagepar Mathieu Lagarde » 15 Sep 2020, 14:45

Bonjour,
Un grand merci à vous deux. J'ai pu adapter le script de Sébastien Rochette avec mon shapefile. C'est parfait ! :)
Bonne journée,
Mat'


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Gabriel Terraz, Google [Bot] et 1 invité