¿Cómo visualizar datos azimutales con incertidumbres?

Estoy tratando de hacer una figura que muestre datos azimutales con un rango diferente de incertidumbres en cada punto. Esta figura de la vieja escuela de un artículo de 1991 captura la idea de la "trama de pajarita" a la que apunto:

¿Alguna sugerencia sobre cómo podría hacer una figura similar? Soy un novato relativo en GIS, pero tengo acceso a ArcGIS a través de mi universidad. Mi experiencia con Arc se ha limitado a hacer mapas geológicos, por lo que no he tenido que hacer nada demasiado exótico.

He examinado las opciones de simbología en Arc y QGIS pero no he visto ninguna configuración que pensé que haría el trabajo. Tenga en cuenta que esto no es solo una cuestión de girar símbolos en forma de corbatín por acimut; El rango angular de cada "corbatín" debe ser diferente.

Yo lo clasificaría como mis conocimientos de Python 'intermedio fuerte' y mis habilidades como R 'intermedio bajo', así que no soy reacio a la piratería algo juntos con matplotliby mpl_toolkits.basemapo bibliotecas similares, si es necesario. Pero pensé en pedir consejo aquí antes de seguir ese camino en caso de que haya una solución más fácil de GIS-land de la que simplemente no haya oído hablar.

Esto requiere un tipo de "cálculo de campo" en el que el valor calculado (basado en la latitud, longitud, acimut central, incertidumbre y distancia) es la forma del corbatín en lugar de un número. Debido a que dichas capacidades de cálculo de campo se hicieron mucho más difíciles en la transición de ArcView 3.x a ArcGIS 8.xy nunca se han restaurado por completo, hoy en día usamos scripts en Python, R o lo que sea: pero el proceso de pensamiento sigue siendo el mismo.

Ilustraré con Rcódigo de trabajo . En su núcleo está el cálculo de una forma de corbatín, que por lo tanto encapsulamos como una función. La función es realmente muy simple: para generar los dos arcos en los extremos del arco, necesita trazar una secuencia a intervalos regulares (de acimut). Esto requiere la capacidad de encontrar las coordenadas (lon, lat) de un punto en función del inicio (lon, lat) y la distancia recorrida. Eso se hace con la subrutina goto, donde ocurre todo el levantamiento aritmético pesado. El resto es solo preparar todo para aplicar gotoy luego aplicarlo.

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

Esto está destinado a aplicarse a una tabla cuyos registros ya debe tener de alguna forma: cada uno de ellos proporciona la ubicación, el acimut, la incertidumbre (como un ángulo a cada lado) y (opcionalmente) una indicación de qué tan grande para hacer el corbata de moño. Simulemos una tabla de este tipo ubicando 1,000 bowties en todo el hemisferio norte:

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

En este punto, las cosas son casi tan simples como cualquier cálculo de campo. Aquí está:

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

(Las pruebas de tiempo indican que Rpuede producir alrededor de 25,000 vértices por segundo. Por defecto, hay un vértice para cada grado de acimut, que el usuario puede configurar mediante el epsargumento bowtie).

Puede hacer un diagrama simple de los resultados en Rsí mismo como un control:

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

Para crear una salida de archivo de forma para importar a un SIG, use el shapefiles paquete:

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

Ahora puede proyectar los resultados, etc. Este ejemplo utiliza una proyección estereográfica del hemisferio norte y las pajaritas están coloreadas por cuantiles de la incertidumbre. (Si observa con mucho cuidado a 180 / -180 grados de longitud, verá dónde este SIG ha recortado las pajaritas que cruzan esta línea. Esa es una falla común con los SIG; no refleja un error en el Rcódigo mismo).