2

I have a data frame of lat/long coordinates and a polygon, which represents a coastline. I am trying to find the distance between each point and the nearest coastline feature. I would like to end up with an output data frame that includes columns for my original lat/long values and a new distance column.

I have attempted to use the gDistance function after reading answers to similar questions online, but I think that I am missing a few steps and am having trouble figuring it out. At present, I only end up with a single distance value. I am quite new to R and would really appreciate any help that anyone might be able to give.

Thanks!

#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: Using AEF's code alterations below (for the for loop step), I am able to get gDistance values for each row, however the output distances are not correct (see below). According to arcGIS they should be between 4-37km, not >500km. Any thoughts on what I am doing wrong here? My land polygon and points are both in the same projection.

gDistance output

  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

The correct distances (calculated in GIS)

  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
PB_Dire
  • 23
  • 5
  • Does [this](http://stackoverflow.com/questions/16448402/distance-of-point-feature-to-nearest-polygon-in-r) or [this](http://stackoverflow.com/questions/15294343/calculating-the-distance-between-polygon-and-point-in-r) help, or did you already look at these? – eipi10 Sep 01 '16 at 18:08
  • @ eipi10 - Thanks for the links. Yes I did read them but they didn't help me to solve the problem. Thanks. – PB_Dire Sep 02 '16 at 11:46

1 Answers1

1

I think you get only one distance value because you overwrite g in every iteration of your for-loop. I do however not not know if this is the only problem because I cannot reproduce your issue without suitable data. Try changing the last loop to this:

g = rep(NA, dim(spdf)[1])
for(i in 1:dim(spdf)[1]){
  g[i] = gDistance(spdf[i,],land_poly)
}
AEF
  • 5,408
  • 1
  • 16
  • 30
  • @ AEF: thanks for the help! The for loop now seems to work, however I still end up with very strange output values (i.e. I should end up with values ~15km and instead get >500km). I have updated my question to include point data, but I cannot upload files to stackoverflow so am not able to provide a shapefile. I tried to use a map from the map package to enable reproducibility but map objects don't work in gDistance. Any further help would be much appreciated. Thanks again! – PB_Dire Sep 02 '16 at 11:45