Bagaimana cara memvisualisasikan data azimut dengan ketidakpastian?

Saya mencoba membuat gambar yang menunjukkan data azimut dengan rentang ketidakpastian yang berbeda di setiap titik. Figur sekolah tua dari makalah tahun 1991 ini menangkap ide "plot bowtie" yang saya tuju:

Ada saran tentang bagaimana saya bisa membuat angka yang sama? Saya seorang pemula di bidang GIS, tetapi saya memiliki akses ke ArcGIS melalui universitas saya. Pengalaman Arc saya terbatas pada membuat peta geologis, jadi saya tidak perlu melakukan sesuatu yang terlalu eksotis.

Saya sudah mencari-cari di opsi simbologi di Arc dan QGIS tetapi belum melihat pengaturan yang saya pikir akan melakukan pekerjaan. Perhatikan bahwa ini bukan hanya masalah memutar simbol berbentuk bowtie oleh azimuth; kisaran sudut setiap "bowtie" harus berbeda.

Saya akan menilai keterampilan Python saya sebagai 'perantara kuat' dan keterampilan R saya sebagai 'perantara menengah', jadi saya tidak segan untuk meretas sesuatu bersama-sama dengan matplotlibdan mpl_toolkits.basemapatau perpustakaan serupa jika perlu. Tapi saya pikir saya akan meminta saran di sini sebelum turun jalan kalau-kalau ada solusi yang lebih mudah dari GIS-tanah yang saya belum pernah dengar.

Ini membutuhkan semacam "perhitungan lapangan" di mana nilai dihitung (berdasarkan pada lintang, bujur, azimuth pusat, ketidakpastian, dan jarak) adalah bentuk bowtie daripada angka. Karena kemampuan perhitungan lapangan seperti itu dibuat jauh lebih sulit dalam transisi dari ArcView 3.x ke ArcGIS 8.x dan belum pernah sepenuhnya pulih, saat ini kami menggunakan scripting dengan Python, R, atau apa pun sebagai gantinya: tetapi proses pemikiran masih merupakan sama.

Saya akan mengilustrasikannya dengan Rkode kerja . Pada intinya adalah perhitungan bentuk bowtie, yang oleh karena itu kami merangkum sebagai fungsi. Fungsi ini benar-benar sangat sederhana: untuk menghasilkan dua busur di ujung haluan, perlu melacak urutan secara berkala (dari azimuth). Ini membutuhkan kemampuan untuk menemukan koordinat (lon, lat) dari suatu titik berdasarkan pada permulaan (lon, lat) dan jarak yang dilalui. Itu dilakukan dengan subrutin goto, di mana semua pengangkatan aritmetika berat terjadi. Sisanya hanya menyiapkan segala sesuatu untuk diterapkan gotodan kemudian menerapkannya.

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

Ini dimaksudkan untuk diterapkan pada tabel yang catatannya harus sudah Anda miliki dalam beberapa bentuk: masing-masing memberikan lokasi, azimuth, ketidakpastian (sebagai sudut ke setiap sisi), dan (opsional) indikasi seberapa besar untuk membuat dasi kupu-kupu. Mari kita simulasikan tabel seperti itu dengan menempatkan 1.000 dasi di seluruh belahan bumi utara:

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

Pada titik ini, semuanya hampir sesederhana perhitungan bidang apa pun. Ini dia:

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

(Tes waktu menunjukkan bahwa Rdapat menghasilkan sekitar 25.000 simpul per detik. Secara default, ada satu simpul untuk setiap derajat azimuth, yang dapat diatur pengguna melalui epsargumen ke bowtie.)

Anda dapat membuat plot sederhana dari hasil Ritu sendiri sebagai cek:

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

Untuk membuat keluaran shapefile untuk diimpor ke GIS, gunakan shapefilespaket:

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

Sekarang Anda dapat memproyeksikan hasilnya, dll. Contoh ini menggunakan proyeksi stereografi belahan bumi utara dan dasi kupu-kupu diwarnai oleh kuantil ketidakpastian. (Jika Anda melihat dengan sangat hati-hati pada 180 / -180 derajat bujur Anda akan melihat di mana GIS ini telah memotong dasi kupu-kupu yang melintasi garis ini. Itu adalah kesalahan umum dengan GIS; itu tidak mencerminkan bug dalam Rkode itu sendiri.)