태그 보관물: r

r

불확실성으로 방위각 데이터를 시각화하는 방법은 무엇입니까? 통해 ArcGIS에 액세스 할 수 있습니다. 나의

각 지점에서 다른 범위의 불확실성을 가진 방위 데이터를 보여주는 그림을 만들려고합니다. 1991 년 논문의이 구식 그림은 내가 목표로하는 “bowtie plot”아이디어를 포착합니다.

Hillhouse and Wells, 1991. "아리조나, 캘리포니아, 네바다에있는 낮은 Miocene 복숭아 스프링 응회암의 자성 직물, 흐름 방향 및 소스 영역"

비슷한 인물을 만드는 방법에 대한 제안 사항이 있습니까? 저는 GIS 관련 신입생이지만 대학을 통해 ArcGIS에 액세스 할 수 있습니다. 나의 Arc 경험은 지리지도를 작성하는 것으로 제한되었으므로 너무 이국적인 일을 할 필요가 없었습니다.

Arc 및 QGIS의 기호 옵션에서 펑크했지만 작업을 수행 할 것으로 생각되는 설정을 보지 못했습니다. 이것은 단지 나비 모양의 기호를 방위각으로 회전시키는 문제가 아니라는 점에 유의하십시오. 각 “bowtie”의 각도 범위는 달라야합니다.

필자는 파이썬 기술을 ‘강한 중간체’로, R 기술을 ‘낮은 중간체’로 평가하므로 필요한 경우 비슷한 라이브러리 matplotlib와 함께 mpl_toolkits.basemap또는 유사한 라이브러리 를 해킹하는 것을 피하지 않습니다 . 그러나 방금 들어 본 적이없는 GIS 랜드의 더 쉬운 솔루션이있는 경우 길을 가기 전에 여기에 대한 조언을 구할 것이라고 생각했습니다.



답변

여기에는 계산 된 값 (위도, 경도, 중심 방위각, 불확실성 및 거리에 기초한 값)이 숫자가 아닌 나비 모양 인 일종의 “필드 계산”이 필요합니다. 이러한 필드 계산 기능은 ArcView 3.x에서 ArcGIS 8.x로 전환하는 데 훨씬 어려워졌으며 완전히 복원 된 적이 없기 때문에 오늘날에는 Python, R 또는 그 밖의 스크립트를 사용합니다. 그러나 사고 과정은 여전히 같은.

작업 R코드로 설명하겠습니다 . 핵심은 나비 넥타이 모양의 계산입니다. 따라서 함수로 캡슐화합니다. 이 기능은 매우 간단합니다. 활의 끝에 두 개의 호를 생성하려면 규칙적인 간격 (방위각)으로 시퀀스를 추적해야합니다. 이를 위해서는 시작 (lon, lat) 및 이송 거리를 기준으로 점의 (lon, lat) 좌표를 찾는 기능이 필요합니다. 이것은 goto모든 무거운 산술 리프팅이 발생 하는 서브 루틴으로 이루어집니다 . 나머지는 적용 할 모든 것을 준비한 goto다음 적용하는 것입니다.

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

이것은 레코드가 이미 어떤 형식으로되어 있어야하는 테이블에 적용하기위한 것입니다. 각 레코드는 위치, 방위각, 불확실성 (각 측면에 대한 각도) 및 (선택적) 나비 넥타이. 북반구에 1,000 개의 나비를 배치하여 이러한 테이블을 시뮬레이션 해 보겠습니다.

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

이 시점에서 모든 것은 필드 계산만큼 간단합니다. 여기있어:

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

(타이밍 테스트는 R초당 약 25,000 개의 정점을 생성 할 수 있음을 나타냅니다 . 기본적으로 각 방위각에 대해 하나의 정점이 있으며,이 eps인수 는를 통해 사용자가 설정할 수 bowtie있습니다.)

R검사 자체로 결과의 간단한 플롯을 만들 수 있습니다 .

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

R로 플롯

GIS로 가져 오기 위해 shapefile 출력을 작성하려면 shapefiles패키지를 사용하십시오 .

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

이제 결과 등을 투영 할 수 있습니다.이 예에서는 북반구의 입체 투영법을 사용하고 나비 넥타이는 불확실성의 수에 의해 채색됩니다. (180 / -180도 경도를주의 깊게 살펴보면이 GIS가이 선을 가로 지르는 나비 넥타이를 어디에서 잘랐는 지 알 수 있습니다. 이는 GIS의 일반적인 결함이며 R코드 자체 의 버그는 반영하지 않습니다 .)

ArcView에서 플롯


답변