Como visualizar dados azimutais com incertezas?

Estou tentando fazer uma figura mostrando dados azimutais com um intervalo diferente de incertezas em cada ponto. Esta figura da velha escola de um artigo de 1991 captura a idéia do "enredo da gravata borboleta" que pretendo:

Alguma sugestão sobre como eu poderia fazer uma figura semelhante? Sou novato em assuntos relacionados ao GIS, mas tenho acesso ao ArcGIS através da minha universidade. Minha experiência com o arco foi limitada à criação de mapas geológicos, por isso não tive que fazer nada muito exótico.

Analisei as opções de simbologia no Arc e no QGIS, mas não vi nenhuma configuração que eu pensasse que faria o trabalho. Observe que não se trata apenas de girar símbolos em forma de gravata borboleta pelo azimute; o alcance angular de cada "gravata borboleta" precisa ser diferente.

Eu classificaria minhas habilidades em Python como 'intermediário forte' e minhas habilidades de R como 'intermediário baixo', por isso não sou averso a invadir algo junto com matplotlibe / mpl_toolkits.basemapou bibliotecas semelhantes, se necessário. Mas pensei em solicitar aconselhamento aqui antes de seguir esse caminho, caso haja uma solução mais fácil da área de SIG que eu simplesmente não tenha ouvido falar.

Isso requer um tipo de "cálculo de campo" no qual o valor calculado (com base em latitude, longitude, azimute central, incerteza e distância) é a forma da gravata borboleta em vez de um número. Como esses recursos de cálculo de campo foram muito mais difíceis na transição do ArcView 3.x para o ArcGIS 8.xe nunca foram totalmente restaurados, hoje em dia usamos scripts em Python, R ou qualquer outra coisa: o processo de pensamento ainda é o mesmo.

Ilustrarei com Rcódigo de trabalho . Em sua essência, está o cálculo de uma forma de gravata borboleta, que, portanto, encapsulamos como uma função. A função é realmente muito simples: para gerar os dois arcos nas extremidades do arco, é necessário traçar uma sequência em intervalos regulares (de azimute). Isso requer a capacidade de encontrar as coordenadas (lon, lat) de um ponto com base na partida (lon, lat) e na distância percorrida. Isso é feito com a sub-rotina goto, onde ocorre todo o levantamento aritmético pesado. O resto é apenas preparar tudo para aplicar gotoe depois aplicá-lo.

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

Isso deve ser aplicado a uma tabela cujos registros você já deve ter de alguma forma: cada um deles fornece a localização, o azimute, a incerteza (como um ângulo para cada lado) e (opcionalmente) uma indicação de quão grande deve ser o tamanho do arquivo. gravata-borboleta. Vamos simular essa tabela situando 1.000 gravatas-borboleta em todo o hemisfério 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)
  ))

Neste ponto, as coisas são quase tão simples quanto qualquer cálculo de campo. Aqui 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)))))

(Os testes de tempo indicam que Rpode produzir cerca de 25.000 vértices por segundo. Por padrão, há um vértice para cada grau de azimute, que pode ser definido pelo usuário pelo epsargumento de bowtie.)

Você pode fazer um gráfico simples dos resultados em Rsi como uma verificação:

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 criar uma saída shapefile para importar para um GIS, use o shapefilespacote:

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

Agora você pode projetar os resultados, etc. Este exemplo usa uma projeção estereográfica do hemisfério norte e os laços são coloridos por quantis da incerteza. (Se você observar com muito cuidado a 180 / -180 graus de longitude, verá onde esse GIS cortou as gravatas-borboleta que cruzam essa linha. Essa é uma falha comum dos GIS; não reflete um erro no Rpróprio código.)