2

何度も使用geosphere::areaPolygon()して成功しましたが、今問題が発生しました。

別のポリゴン sp1 を含むポリゴン sp2 があります。したがって、sp2 は sp1 よりも大きな領域を持つ必要があります。で各面積を計算するとareaPolygon()、逆の結果が得られます。

areasp1 = 10133977
areasp2 = 9858811

gSymdifference を使用して対称的な異なるポリゴン sp3 を見つけました。

areasp3 = 275165.4

sp1 : 赤の破線
sp2 : 黒の線
sp3 : 緑の点線 画像はこちら

私が使用したコード:

library(sp)
library(geosphere)
library(rgeos)
library(raster)

s1 <- data.frame(lon= c(130.4327, 129.85127, 121.00775, 84.9176, 60.40123, 
58.97929, 58.55841, 94.95237),lat= c(30.25074, 29.83075, 23.7992, 28.25964, 
36.89905, 37.72305, 52.58793, 43.47459))

s2 <- data.frame(lon= c(130.4327, 129.85127, 121.00775, 54.35377, 58.55841, 
94.95237),lat= c(30.25074, 29.83075, 23.7992, 31.61798, 52.58793, 43.47459))

sp1 <- SpatialPolygons(list(Polygons(list(Polygon(s1)),1)))
crs(sp1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

sp2 <- SpatialPolygons(list(Polygons(list(Polygon(s2)),1)))
crs(sp2) <-CRS ("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

areasp1 <- areaPolygon(sp1)/10^6    # to get area in square km
areasp2 <- areaPolygon(sp2)/10^6

sp3 <- gSymdifference(sp1,sp2,checkvalidity = TRUE)
areasp3 <- areaPolygon(sp3)/10^6 

すべてのポリゴンは、 の自己交差またはその他の問題についてテストされているため、ここでgIsValid()言及されている問題とは関係ありません。

sp1がsp2よりも面積が大きい理由について何か考えはありますか?

注:areasp2 +を合計areasp3すると、ほぼ に等しくなりareasp1ます。これが偶然かどうかはわかりません。

4

1 に答える 1