Jak wizualizować dane azymutalne z niepewnością?

Próbuję zrobić rysunek przedstawiający dane azymutalne z innym zakresem niepewności w każdym punkcie. Ta oldschoolowa postać z gazety z 1991 roku oddaje ideę „fabuły muszki”, do której dążę:

Wszelkie sugestie, w jaki sposób mogę zrobić podobną liczbę? Jestem względnym nowicjuszem w GIS, ale mam dostęp do ArcGIS przez mój uniwersytet. Moje doświadczenie Arc jest ograniczone do tworzenia map geologicznych, więc nie musiałem robić nic zbyt egzotycznego.

Przeszukiwałem opcje symboli w Arc i QGIS, ale nie widziałem żadnych ustawień, które według mnie byłyby odpowiednie. Zauważ, że nie chodzi tylko o obracanie symboli w kształcie muszki o azymut; zakres kątowy każdego „muszki” musi być inny.

Chciałbym ocenić moje umiejętności Python jako „silne pośredni” i moich umiejętności R jak „niska pośrednia”, więc nie jestem przeciwny hacking coś razem z matplotlibi mpl_toolkits.basemaplub podobnych bibliotek w razie potrzeby. Ale pomyślałem, że zasięgnę porady tutaj, zanim pójdę tą drogą, na wypadek, gdyby było łatwiejsze rozwiązanie z GIS-landu, o którym po prostu nie słyszałem.

Wymaga to pewnego rodzaju „obliczenia pola”, w którym obliczona wartość (oparta na szerokości i długości geograficznej, azymucie centralnym, niepewności i odległości) jest kształtem muszki, a nie liczbą. Ponieważ takie możliwości obliczania pól były znacznie utrudnione przy przejściu z ArcView 3.x do ArcGIS 8.x i nigdy nie zostały w pełni przywrócone, obecnie używamy skryptów w Pythonie, R itp., Ale proces myślowy jest nadal podobnie.

Zilustruję działającym Rkodem. U ich podstaw leży obliczenie kształtu muszki, który dlatego enkapsulujemy jako funkcję. Ta funkcja jest naprawdę bardzo prosta: aby wygenerować dwa łuki na końcach łuku, należy wyśledzić sekwencję w regularnych odstępach czasu (azymutu). Wymaga to zdolności do znalezienia współrzędnych (lon, lat) punktu na podstawie początkowej (lon, lat) i pokonanej odległości. Odbywa się to z podprogramem goto, w którym odbywa się wszystkie ciężkie arytmetyczne podnoszenie. Reszta przygotowuje wszystko do zastosowania, gotoa następnie stosuje.

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

Ma to zostać zastosowane do tabeli, której rekordy musisz już mieć w jakiejś formie: każda z nich podaje położenie, azymut, niepewność (jako kąt z każdej strony) i (opcjonalnie) wskazanie, jak duże jest wykonanie muszka. Symulujmy taki stół, umieszczając 1000 muszek na całej półkuli północnej:

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

W tym momencie rzeczy są prawie tak proste, jak każde obliczenie pola. Oto on:

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

(Testy czasowe wskazują, że Rmoże wytwarzać około 25 000 wierzchołków na sekundę. Domyślnie dla każdego stopnia azymutu istnieje jeden wierzchołek, który można ustawić za pomocą epsargumentu do bowtie.)

Możesz wykonać prosty wykres wyników Rsam w sobie jako kontrolę:

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

Aby utworzyć plik wyjściowy shapefile do importowania do GIS, użyj shapefilespakietu:

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

Teraz możesz wyświetlać wyniki itp. W tym przykładzie wykorzystano stereograficzną projekcję półkuli północnej, a muszki są zabarwione kwantylami niepewności. (Jeśli przyjrzysz się bardzo uważnie długości geograficznej 180 / -180 stopni, zobaczysz, gdzie GIS przeciął muszki, które przecinają tę linię. Jest to powszechna wada w przypadku GISes; nie odzwierciedla błędu w samym Rkodzie.)