Wie kann man azimutale Daten mit Unsicherheiten visualisieren?

Ich versuche, eine Zahl zu erstellen, die azimutale Daten mit unterschiedlichen Unsicherheiten an jedem Punkt zeigt. Diese Oldschool-Figur aus einem Papier von 1991 fängt die Idee der "Fliege-Handlung" ein, die ich anstrebe:

Irgendwelche Vorschläge, wie ich eine ähnliche Figur machen könnte? Ich bin ein relativer Neuling bei GIS, aber ich habe über meine Universität Zugriff auf ArcGIS. Meine Arc-Erfahrung beschränkte sich darauf, geologische Karten zu erstellen, daher musste ich nichts zu Exotisches tun.

Ich habe mich in den Symbologieoptionen in Arc und QGIS umgesehen, aber keine Einstellungen gesehen, von denen ich dachte, dass sie den Job erledigen würden. Beachten Sie, dass es nicht nur darum geht, bogenförmige Symbole durch Azimut zu drehen. Der Winkelbereich jeder "Fliege" muss unterschiedlich sein.

Ich würde meine Python-Fähigkeiten als "starke Mittelstufe" und meine R-Fähigkeiten als "niedrige Mittelstufe" bewerten, daher bin ich nicht abgeneigt, bei Bedarf etwas zusammen mit matplotlibund mpl_toolkits.basemapoder ähnlichen Bibliotheken zu hacken . Aber ich dachte, ich würde hier Rat einholen, bevor ich diesen Weg gehe, falls es eine einfachere Lösung aus dem GIS-Land gibt, von der ich gerade noch nichts gehört habe.

Dies erfordert eine Art "Feldberechnung", bei der der berechnete Wert (basierend auf Breite, Länge, zentralem Azimut, Unsicherheit und Entfernung) eher die Fliege als eine Zahl ist. Da solche Feldberechnungsfunktionen beim Übergang von ArcView 3.x zu ArcGIS 8.x erheblich erschwert wurden und nie vollständig wiederhergestellt wurden, verwenden wir heutzutage Skripte in Python, R oder was auch immer. Der Denkprozess ist jedoch immer noch der gleich.

Ich werde mit Arbeitscode veranschaulichen R. Im Zentrum steht die Berechnung einer Fliege, die wir daher als Funktion einkapseln. Die Funktion ist wirklich sehr einfach: Um die beiden Bögen an den Enden des Bogens zu erzeugen, muss in regelmäßigen Abständen (von Azimut) eine Sequenz nachgezeichnet werden. Dies erfordert die Fähigkeit, die (lon, lat) -Koordinaten eines Punktes basierend auf dem Start (lon, lat) und der zurückgelegten Strecke zu finden. Dies geschieht mit dem Unterprogramm goto, bei dem das gesamte schwere arithmetische Heben erfolgt. Der Rest bereitet einfach alles für die Bewerbung vor gotound wendet es dann an.

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

Dies soll auf eine Tabelle angewendet werden, deren Datensätze Sie bereits in irgendeiner Form haben müssen: Jede von ihnen gibt die Position, den Azimut, die Unsicherheit (als Winkel zu jeder Seite) und (optional) einen Hinweis darauf an, wie groß die sein sollen Krawatte. Simulieren wir eine solche Tabelle, indem wir 1.000 Bowties auf der gesamten Nordhalbkugel platzieren:

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

Zu diesem Zeitpunkt sind die Dinge fast so einfach wie bei jeder Feldberechnung. Hier ist es:

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

(Timing-Tests zeigen, dass Rungefähr 25.000 Scheitelpunkte pro Sekunde erzeugt werden können. Standardmäßig gibt es einen Scheitelpunkt für jeden Azimutgrad, der vom Benutzer über das epsArgument auf einstellbar ist bowtie.)

Sie können Rzur Überprüfung eine einfache Darstellung der Ergebnisse selbst erstellen:

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

Verwenden Sie das shapefilesPaket, um eine Shapefile-Ausgabe für den Import in ein GIS zu erstellen:

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

Jetzt können Sie die Ergebnisse usw. projizieren. In diesem Beispiel wird eine stereografische Projektion der nördlichen Hemisphäre verwendet, und die Fliegen werden durch Quantile der Unsicherheit gefärbt. (Wenn Sie sich den Längengrad von 180 / -180 Grad genau ansehen, werden Sie sehen, wo dieses GIS die Fliegen abgeschnitten hat, die diese Linie kreuzen. Dies ist ein häufiger Fehler bei GISes; es spiegelt keinen Fehler im RCode selbst wider .)