[Résolu] Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

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

jean lobry
Messages : 733
Enregistré le : 17 Jan 2008, 20:00
Contact :

[Résolu] Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

Messagepar jean lobry » 24 Oct 2018, 19:21

Bonjour à tous,

je suis tout béotien pour les statistiques spatiales. J'essaye de reproduire la figure page 9 de ce rapport à propos des bébés sans main dans le département de l'Ain. Voici où j'en suis avec le code reproductible suivant :

Code : Tout sélectionner

load(url("http://pbil.univ-lyon1.fr/R/donnees/qrbb/AinCol.Rda"))
par(mar = c(0, 0, 2, 0) + 0.1)
plot(ain, col = ain$col)
grid(lty = 1)
main <- "Localisation des 7 cas d'agénésie des membres supérieurs\ndans le département de l'Ain entre 2009 et 2014, inclus"
title(main = main, line = -1)
par(lend = "butt")
legend("bottomleft", inset = 0, lwd = 10,
       legend = c("Cas d'agénésie", "Cluster d'agénésie de la main"),
       col = c("yellow", "lightblue"))

Je ne suis pas satisfait du rendu, j'aimerais avoir l'équivalent de asp = 1 pour les plots standards: mon grid() devrait me faire des carrés en première approximation. Je me doute bien que c'est compliqué dès lors que l'on projette un ellipsoïde sur un plan. Est-ce que quelqu'un aurait une solution simple pour moi ?

Amicalement,

jean

François Bonnot
Messages : 537
Enregistré le : 10 Nov 2004, 15:19
Contact :

Re: Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

Messagepar François Bonnot » 25 Oct 2018, 08:20

Bonjour Jean,
je suis tout béotien pour les statistiques spatiales

Nous est-il permis d'en douter ?

j'aimerais avoir l'équivalent de asp = 1

Le paramètre asp fonctionne (il suffit d'essayer asp=2 pour s'en persuader) mais ce n'est pas suffisant parce que la fonction plot ne fait qu'utiliser les coordonnées stockées dans l'objet.
Il n'y a pas de système de coordonnées affecté à l'objet ain:

Code : Tout sélectionner

> proj4string(ain)
[1] NA

En faisant l'hypothèse que les coordonnées sont en WGS84, la ligne suivante suivie de plot() donne une image qui ressemble à celle du rapport:

Code : Tout sélectionner

proj4string(ain) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
François

François Bonnot
Messages : 537
Enregistré le : 10 Nov 2004, 15:19
Contact :

Re: Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

Messagepar François Bonnot » 25 Oct 2018, 08:46

J'ajouterai : plutôt que grid, on peut utiliser le code suivant pour tracer les vrais méridiens et parallèles (inspiré je crois d'un code de Renaud).

Code : Tout sélectionner

library(rgdal)
projstring <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
meri1 <- -10; meri2 <- 10; para1 <- 40; para2 <- 60
meri1 <- -10; meri2 <- 15; para1 <- 40; para2 <- 60; step <- 0.1
## méridiens
meri <- seq(meri1, meri2, step)
y <- seq(para1, para2, length = 100)
for(i in seq(length(meri))){
    x <- rep(meri[i], 100)
    lines(project(cbind(x, y), projstring), col = ifelse(x==0, "red", "grey"), lwd = .5)
}

## parallèles
para <- seq(para1, para2, step)
x <- seq(meri1, meri2, length = 100)
for(i in seq(length(para))){
    y <- rep(para[i], 100)
    lines(project(cbind(x, y), projstring), col = "grey", lwd = .5)
}
box()
François

jean lobry
Messages : 733
Enregistré le : 17 Jan 2008, 20:00
Contact :

Re: Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

Messagepar jean lobry » 25 Oct 2018, 16:55

Bonjour François,
François Bonnot a écrit :Nous est-il permis d'en douter ?

Non, non, les seules statistiques spatiales que j'ai jamais pratiquées, c'est juste du 1D, genre les chirochores le long des génomes bactériens, jamais des coordonnées GPS de la vraie vie. Mais j'ai bien envie de m'y mettre, c'est quand même sympa que de produire des cartes.

Merci beaucoup pour toutes tes pistes, j'explore tout ça et je reviens vers le forum.

Bien amicalement,

jean

jean lobry
Messages : 733
Enregistré le : 17 Jan 2008, 20:00
Contact :

Re: [Résolu] Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

Messagepar jean lobry » 26 Oct 2018, 17:38

Bonjour à tous,

ci-après le code reproductible et commenté qui explique pourquoi ma question était idiote. Merci à François pour m'avoir mis sur la piste !

Bien amicalement,

jean

Code : Tout sélectionner

load(url("http://pbil.univ-lyon1.fr/R/donnees/qrbb/AinCol.Rda"))
library("sp")
proj4string(ain) # NA, pas glop !
#
# Je suis parti du fichier "communes-20180101.shp" donnant le découpage
# administratif communal français issu d'OpenStreetMap dont j'ai extrait les
# communes du département de l'Ain. Les coordonnées sont en WGS84.
# Le premier truc à faire c'est de le préciser
#
CRS1 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
proj4string(ain) <- CRS1
#
# La clef pour comprendre ce qu'il se passe c'est de tracer les vrais
# méridiens et parallèles. Je bricole une fonction kleenex pour ça.
#
addlonlat <- function(xy, CRS1, CRS2){
  lon.min <- floor(min(xy[ , 1]))
  lon.max <- ceiling(max(xy[ , 1]))
  lon.step <- 0.1
  lon.seq <- seq(from = lon.min, to = lon.max, by = lon.step)
  lon.col <- rep("grey", length(lon.seq))
  lon.col[which(sapply(lon.seq - floor(lon.seq), all.equal, 0) == "TRUE")] <- "red"
  lat.min <- floor(min(xy[ , 2]))
  lat.max <- ceiling(max(xy[ , 2]))
  lat.step <- 0.1
  lat.seq <- seq(from = lat.min, to = lat.max, by = lat.step)
  lat.col <- rep("grey", length(lat.seq))
  lat.col[which(sapply(lat.seq - floor(lat.seq), all.equal, 0) == "TRUE")] <- "red"
 
  for(i in 1:length(lon.seq))
  {
    myLine <- Line(coords = cbind(c(lon.seq[i], lon.seq[i]), c(lat.min, lat.max)))
    myLines <- Lines(slinelist = list(myLine), ID = "ID")
    mySpatialLines <- SpatialLines(list(myLines), proj4string = CRS1)
    mySpatialLines <- spTransform(mySpatialLines, CRS2)
    plot(mySpatialLines, add = TRUE, col = lon.col[i])
  }
 
  for(j in 1:length(lat.seq))
  {
    myLine <- Line(coords = cbind(c(lon.min, lon.max), c(lat.seq[j], lat.seq[j])))
    myLines <- Lines(slinelist = list(myLine), ID = "ID")
    mySpatialLines <- SpatialLines(list(myLines), proj4string = CRS1)
    mySpatialLines <- spTransform(mySpatialLines, CRS2)
    plot(mySpatialLines, add = TRUE, col = lat.col[j])
  }
}
#
# Je teste ma fonction kleenex avec une projection bien tordue pour vérifier
# que les longitudes et latitudes sont projetées correctement.
#
xy <- coordinates(ain) # les coordonnées de départ
CRS2 <- CRS("+proj=gnom +lat_0=90 +lon_0=-60 +no_def")
ain.gnom <- spTransform(ain, CRS2)
par(mar = rep(0, 4) + 0.1)
plot(ain.gnom, col = ain.gnom$col, asp = 1)
addlonlat(xy, CRS1, CRS2)
#
# Je reviens à la projection de base, c'est un peu étiré dans le sens Nord-Sud,
# comme dans la carte du rapport
#
plot(ain, col = ain$col)
addlonlat(xy, CRS1, CRS1)
#
# Mais si on précise asp = 1, tout va bien
#
plot(ain, col = ain$col, asp = 1)
addlonlat(xy, CRS1, CRS1)
#
# Conclusion, asp = 1 marche très bien. Ne pas se fier à grid()
#

François Bonnot
Messages : 537
Enregistré le : 10 Nov 2004, 15:19
Contact :

Re: [Résolu] Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

Messagepar François Bonnot » 29 Oct 2018, 08:45

Mais si on précise asp = 1, tout va bien

Jean,
Je n'en suis pas si certain.
Je me demande si le paramètre asp=1 n'est pas en trop, en effet il a pour effet de rendre la grille carrée, i.e. la longueur d'un arc de méridien égale à celle d'un arc de parallèle de même angle, au lieu d'un rapport d'environ 0.7 sous une latitude de 45°. Sans asp=1, le rapport de 0.7 semble respecté, ce qui signifierait que plot appliqué à un objet de classe SpatialPolygonsDataFrame gère automatiquement le paramètre asp.
Mon interprétation (les spécialistes rectifieront si nécessaire) est que l'objet ain.gnom, issu d'une vraie projection (gnonomique), a des coordonnées métriques donc représentables avec asp=1 (d'ailleurs la suppression de ce paramètre ne change rien), alors que ain a des coordonnées angulaires qui font que le paramètre asp=1 donne la même longueur au même angle de longitude et de latitude.
Merci pour la fonction addlonlat que chacun pourra réutiliser.
François

jean lobry
Messages : 733
Enregistré le : 17 Jan 2008, 20:00
Contact :

Re: [Résolu] Quel est l'équivalent de asp = 1 pour des SpatialPolygonsDataFrame ?

Messagepar jean lobry » 29 Oct 2018, 14:42

François Bonnot a écrit :
Mais si on précise asp = 1, tout va bien

Je n'en suis pas si certain.
Je me demande si le paramètre asp=1 n'est pas en trop, en effet il a pour effet de rendre la grille carrée, i.e. la longueur d'un arc de méridien égale à celle d'un arc de parallèle de même angle, au lieu d'un rapport d'environ 0.7 sous une latitude de 45°. Sans asp=1, le rapport de 0.7 semble respecté, ce qui signifierait que plot appliqué à un objet de classe SpatialPolygonsDataFrame gère automatiquement le paramètre asp.


Tu as complètement raison, je raisonne faussement comme si la dérive des continents avait ramené le département de l'Ain au niveau de l'équateur terrestre, ça risque de prendre encore un petit bout de temps :-)

Amicalement

jean


Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité