여기에는 계산 된 값 (위도, 경도, 중심 방위각, 불확실성 및 거리에 기초한 값)이 숫자가 아닌 나비 모양 인 일종의 “필드 계산”이 필요합니다. 이러한 필드 계산 기능은 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))

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

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

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

답변