2011-12-20 63 views
1

我有一個包含50多個不同多邊形形狀(代表50多個不同區域)的shapefile和應該存在於其中一個區域的10,000多個數據點。事情是,10,000+個點已經編碼了一個他們應該在的區域,我想知道他們距離這個編碼區域在地理空間距離上有多遠。如何從多邊形計算地理距離?

我目前的方法(代碼如下),它涉及到從sp庫轉換形狀文件到owin對象和使用distfun讓我在lat,long euclidean空間的距離。但我想獲得地理空間距離(最終轉換爲km)。我應該去哪裏?

#basically cribbed from http://cran.r-project.org/web/packages/spatstat/vignettes/shapefiles.pdf (page 9) 
shp <- readShapeSpatial("myShapeFile.shp", proj4string=CRS("+proj=longlat +datum=WGS84")) 
regions <- lapply(slot(shp, "polygons"), function(x) SpatialPolygons(list(x))) 
windows <- lapply(regions, as.owin) 

# need to convert this to geo distance 
distance_from_region <- function(regionData, regionName) { 
    w <- windows[[regionName]] 
    regionData$dists <- distfun(w)(regionData$lat, regionData$long) 
    regionData 
} 
+0

在sp中嘗試spDistsN1,使用longlat = TRUE – mdsumner 2011-12-20 20:11:49

+0

@mdsummer:spDistsN1計算兩點之間的距離。我需要從多邊形到點的距離。我不明白如何做轉換... – prabhasp 2011-12-20 21:47:46

+0

多邊形的座標?質心或某個選定的邊界點(最近的一個?),還是所有的邊界點?你的問題不清楚。 spDists/spDistsN1在沒有選擇地圖投影的情況下提供了您需要的多組點之間的大圓距離(不在兩點之間),所以一旦提取了正確的座標或摘要,我最終會使用它 - 您能否澄清這個問題? – mdsumner 2011-12-20 22:00:39

回答

3

我項目的數據歐幾里得(或接近歐幾里德)座標系統 - 除非你是橫跨一大塊地球的,那麼這是可行的。使用來自maptools或sp或rgdal的spTransform(我忘記了它)並將其轉換爲數據附近的UTM區域。

你也可能與包裝rgeos和gDistance功能做的更好:

gDistance by default returns the cartesian minimum distance 
between the two geometries in the units of the current projection. 

如果你的數據是在一大塊全球然後... ...棘手...... 42

Barry