2016-09-01 14 views
2

У меня есть кадр данных координат lat/long и многоугольник, представляющий собой береговую линию. Я пытаюсь найти расстояние между каждой точкой и ближайшей функцией береговой линии. Я хотел бы получить рамку выходных данных, которая включает столбцы для моих первоначальных значений lat/long и новый столбец расстояния.Поиск ближайших расстояний между вектором точек и многоугольником в R

Я попытался использовать функцию gDistance после прочтения ответов на подобные вопросы в Интернете, но я думаю, что у меня не хватает нескольких шагов, и мне трудно понять это. В настоящее время я получаю только одно значение расстояния. Я совершенно новичок в R и буду очень признателен за любую помощь, которую любой может дать.

Спасибо!

#Load data 
    Locs = structure(list(id = 1:5, Lat = c(29.59679167, 29.43586667, 29.37642222,29.52786111, 30.10603611), Long = c(-81.02547778, -80.92573889, 
-80.97714167, -81.08721667, -80.94368611)), .Names = c("id","Lat", "Long"), class = "data.frame", row.names = c(NA, -5L)) 

#Extract lat/long coordinates 
    xy = Locs[,c("Lat","Long")] 
#Create SpatialPointsDataFrame from xy data and change projection to metres 
    spdf = SpatialPointsDataFrame(coords=xy, data=xy, proj4string = CRS("+proj=aea +zone=17 ellps=WGS84")) 

#Read in shapefile as a spatialdataframe object 
    coast = readOGR(dsn="land data", layer="coast")  
#Transform to AEA (m) projection to match projection of points 
    land_poly = spTransform(coast, CRS("+proj=aea +zone=17 ellps=WGS84")) 

#OR load map from map package (but unfortunately map objects do not work in gDistance) 
    library(maps) 
    library(mapdata) 
    coast2 = map('usa', col = "grey90", fill=TRUE) 

#Calculate distance between each point and the nearest land feature 
    for(i in 1:dim(spdf)[1]){ 
    g = gDistance(spdf[i,],land_poly) 
    } 

EDIT: Использование кода изменения АЭФ в ниже (для для шага цикла), я могу получить значения gDistance для каждой строки, однако выходные расстояния не являются правильными (см. Ниже) Согласно arcGIS они должны быть между 4-37 км, а не> 500 км. Любые мысли о том, что я делаю неправильно здесь? Мой земельный многоугольник и точки находятся в одной проекции.

выход gDistance

id  Lat  Long dist_gDist 
1 1 29.59679 -81.02548 516299.0 
2 2 29.43587 -80.92574 516298.8 
3 3 29.37642 -80.97714 516298.9 
4 4 29.52786 -81.08722 516299.0 
5 5 30.10604 -80.94369 516299.0 

правильные расстояния (рассчитанные в ГИС)

id  Lat  Long dist_arc 
1 1 29.59679 -81.02548 13.630 
2 2 29.43587 -80.92574 15.039 
3 3 29.37642 -80.97714 8.111 
4 4 29.52786 -81.08722 4.784 
5 5 30.10604 -80.94369 36.855 
+0

ли [это] (http://stackoverflow.com/questions/16448402/distance-of-point-feature-to-nearest-polygon-in-r) или [это] (http://stackoverflow.com/questions/15294343/calculating-the-distance-between-polygon-and-point-in-r), или вы уже посмотрели на них? – eipi10

+0

@ eipi10 - Спасибо за ссылки. Да, я прочитал их, но они не помогли мне решить проблему. Благодарю. –

ответ

1

Я думаю, вы получите только одно значение расстояния, потому что вы перезаписать g в каждой итерации вашей for -loop. Однако я не знаю, является ли это единственной проблемой, потому что я не могу воспроизвести вашу проблему без подходящих данных. Попробуйте изменить последнюю петлю на это:

g = rep(NA, dim(spdf)[1]) 
for(i in 1:dim(spdf)[1]){ 
    g[i] = gDistance(spdf[i,],land_poly) 
} 
+0

@ AEF: спасибо за помощь! Кажется, что цикл for работает, но я все еще получаю очень странные выходные значения (т. Е. Я должен получить значения ~ 15 км и вместо этого получить> 500 км). Я обновил свой вопрос, чтобы включить данные точки, но я не могу загружать файлы в stackoverflow, поэтому не могу предоставить шейп-файл. Я попытался использовать карту из пакета карты, чтобы обеспечить воспроизводимость, но объекты карты не работают в gDistance. Любая дополнительная помощь будет высоко оценена. Еще раз спасибо! –

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

  • Нет связанных вопросов^_^