Comment visualiser les données azimutales avec des incertitudes?

J'essaie de faire une figure montrant des données azimutales avec une gamme différente d'incertitudes à chaque point. Cette figure oldschool d'un papier de 1991 capture l'idée de "noeud papillon" que je vise:

Des suggestions sur la façon dont je pourrais faire une figure similaire? Je suis relativement novice dans le domaine des SIG, mais j'ai accès à ArcGIS via mon université. Mon expérience Arc s'est limitée à faire des cartes géologiques, donc je n'ai rien eu de trop exotique.

J'ai fouillé dans les options de symbologie dans Arc et QGIS, mais je n'ai vu aucun paramètre qui, selon moi, ferait l'affaire. Notez que ce n'est pas seulement une question de rotation des symboles en forme de nœud papillon par azimut; la plage angulaire de chaque "noeud papillon" doit être différente.

Je classerais mes compétences en Python comme «intermédiaire fort» et mes compétences R comme «intermédiaire faible», donc je ne suis pas opposé à pirater quelque chose avec matplotlibet / mpl_toolkits.basemapou des bibliothèques similaires si nécessaire. Mais je pensais que je demanderais des conseils ici avant de s'engager dans cette voie au cas où il y aurait une solution plus facile de SIG-terre dont je n'ai tout simplement pas entendu parler.

Cela nécessite une sorte de "calcul sur le terrain" dans lequel la valeur calculée (basée sur une latitude, une longitude, un azimut central, une incertitude et une distance) est la forme du nœud papillon plutôt qu'un nombre. Étant donné que ces capacités de calcul de champ ont été rendues beaucoup plus difficiles lors de la transition d'ArcView 3.x à ArcGIS 8.x et n'ont jamais été entièrement restaurées, nous utilisons aujourd'hui des scripts en Python, R ou autre chose à la place: mais le processus de réflexion est toujours le même.

Je vais illustrer avec du Rcode de travail . Au cœur se trouve le calcul d'une forme de nœud papillon, que nous encapsulons donc en fonction. La fonction est vraiment très simple: pour générer les deux arcs aux extrémités de l'arc, il faut tracer une séquence à intervalles réguliers (d'azimut). Cela nécessite une capacité à trouver les coordonnées (lon, lat) d'un point en fonction du départ (lon, lat) et de la distance parcourue. Cela se fait avec le sous-programme goto, où tout le levage arithmétique lourd se produit. Le reste prépare simplement tout pour postuler goto, puis l'applique.

bowtie <- function(azimuth, delta, origin=c(0,0), radius=1, eps=1) {
  #
  # On entry:
  #   azimuth and delta are in degrees.
  #   azimuth is east of north; delta should be positive.
  #   origin is (lon, lat) in degrees.
  #   radius is in meters.
  #   eps is in degrees: it is the angular spacing between vertices.
  #
  # On exit:
  #   returns an n by 2 array of (lon, lat) coordinates describing a "bowtie" shape.
  #
  # NB: we work in radians throughout, making conversions from and to degrees at the
  #   entry and exit.
  #--------------------------------------------------------------------------------#
  if (eps <= 0) stop("eps must be positive")
  if (delta <= 0) stop ("delta must be positive")
  if (delta > 90) stop ("delta must be between 0 and 90")
  if (delta >= eps * 10^4) stop("eps is too small compared to delta")
  if (origin[2] > 90 || origin[2] < -90) stop("origin must be in lon-lat")
  a <- azimuth * pi/180; da <- delta * pi/180; de <- eps * pi/180 
  start <- origin * pi/180
  #
  # Precompute values for `goto`.
  #
  lon <- start[1]; lat <- start[2]
  lat.c <- cos(lat); lat.s <- sin(lat)
  radius.radians <- radius/6366710
  radius.c <- cos(radius.radians); radius.s <- sin(radius.radians) * lat.c
  #
  # Find the point at a distance of `radius` from the origin at a bearing of theta.
  # http://williams.best.vwh.net/avform.htm#Math
  #
  goto <- function(theta) {
    lat1 <- asin(lat1.s <- lat.s * radius.c + radius.s * cos(theta))
    dlon <- atan2(-sin(theta) * radius.s, radius.c - lat.s * lat1.s)
    lon1 <- lon - dlon + pi %% (2*pi) - pi
    c(lon1, lat1)
  }
  #
  # Compute the perimeter vertices.
  #
  n.vertices <- ceiling(2*da/de)
  bearings <- seq(from=a-da, to=a+da, length.out=n.vertices)
  t(cbind(start,
        sapply(bearings, goto),
          start,
        sapply(rev(bearings+pi), goto),
          start) * 180/pi)
}

Ceci est destiné à être appliqué à une table dont vous devez déjà avoir les enregistrements sous une forme quelconque: chacun d'eux donne l'emplacement, l'azimut, l'incertitude (sous la forme d'un angle de chaque côté) et (éventuellement) une indication de la taille de la nœud papillon. Simulons une telle table en situant 1000 bowties dans tout l'hémisphère nord:

n <- 1000
input <- data.frame(cbind(
  id = 1:n, 
  lon = runif(n, -180, 180),
  lat = asin(runif(n)) * 180/pi,
  azimuth = runif(n, 0, 360),
  delta = 90 * rbeta(n, 20, 70),
  radius = 10^7/90 * rgamma(n, 10, scale=2/10)
  ))

À ce stade, les choses sont presque aussi simples que n'importe quel calcul de champ. C'est ici:

  shapes <- as.data.frame(do.call(rbind, 
         by(input, input$id, 
            function(d) cbind(d$id, bowtie(d$azimuth, d$delta, c(d$lon, d$lat), d$radius, 1)))))

(Les tests de synchronisation indiquent qu'il Rpeut produire environ 25 000 sommets par seconde. Par défaut, il existe un sommet pour chaque degré d'azimut, qui peut être défini par l'utilisateur via l' epsargument to bowtie.)

Vous pouvez faire un simple tracé des résultats en Rsoi comme une vérification:

colnames(shapes) <- c("id", "x", "y")
plot(shapes$x, shapes$y, type="n", xlab="Longitude", ylab="Latitude", main="Bowties")
temp <- by(shapes, shapes$id, function(d) lines(d$x, d$y, type="l", lwd=2, col=d$id))

Pour créer une sortie de fichier de formes à importer dans un SIG, utilisez le shapefilespackage:

require(shapefiles)
write.shapefile(convert.to.shapefile(shapes, input, "id", 5), "f:/temp/bowties", arcgis=T)

Vous pouvez maintenant projeter les résultats, etc. Cet exemple utilise une projection stéréographique de l'hémisphère nord et les nœuds papillons sont colorés par les quantiles de l'incertitude. (Si vous regardez très attentivement à 180 / -180 degrés de longitude, vous verrez où ce SIG a coupé les nœuds papillons qui traversent cette ligne. C'est un défaut commun avec les SIG; il ne reflète pas un bug dans le Rcode lui-même.)