2011-12-20 4 views
1

У меня есть шейп-файл с 50 + разными многоугольными формами (представляющими 50+ разных регионов) и 10 000+ точек данных, которые должны присутствовать в одном из регионов. Дело в том, что 10 000 + баллов уже закодированы с областью, в которой они должны находиться, и я хочу выяснить, насколько они далеко от этой кодированной области на геопространственном расстоянии.Как рассчитать гео-расстояние от многоугольника?

Мой текущий подход (код ниже), который включает в себя преобразование шейпфайлов в owin объектов из sp библиотеки и использования distfun получает меня расстояние в широчайшем, долго евклидове пространства. Но я хотел бы получить геопространственные расстояния (в конечном счете, чтобы перейти в км). Куда мне идти дальше?

#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

Попробуйте spDistsN1 в sp, with longlat = TRUE – mdsumner

+0

@mdsummer: spDistsN1 вычисляет расстояния между двумя точками. Мне нужно расстояние от полигона до точки. Мне сразу не ясно, как сделать преобразование ... – prabhasp

+0

, из которого координаты многоугольника? Центроид или какая-то выбранная граничная точка (ближайшая?) Или все граничные точки? На этот вопрос неясно. spDists/spDistsN1 предоставляет большие расстояния на расстоянии между наборами точек (не только между двумя точками), не выбирая проекцию карты, поэтому в конечном счете я буду использовать, когда будут извлечены правильные координаты или сводка, - можете ли вы уточнить вопрос об этом ? – mdsumner

ответ

3

Я бы проецировать данные евклидову (или вблизи евклидовой) системы координат - если вы не охватывает значительную часть земного шара, то это возможно. Используйте spTransform из maptools или sp или rgdal (я забыл, что) и конвертируем в зону UTM рядом с вашими данными.

Вы также мог бы сделать лучше с пакетом rgeos и функцией gDistance:

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

Если данные на большой кусок земного шара, то ... сложно ... 42 ...

Barry

Смежные вопросы