Search code examples
rgeospatialpolygongeosphere

Geosphere package distance from point to polygon


I'm trying to use the geosphere package in R to get the distance to a polygon from a set of points that lie outside of that polygon.

The polygon is a shapefile of the Antarctic coastline, found here: https://data.bas.ac.uk/items/e6cf8946-e493-4c36-b4f5-58f7a2ee2a74/ and the points are animal tracking data.

I have tried using the syntax specified in the geosphere documentation (https://www.rdocumentation.org/packages/geosphere/versions/1.5-14/topics/dist2Line) which is as follows:

dist2Line(p, line, distfun=distGeo)
#my attempt so far: 
#libraries 
library(rgdal)
library(sf)
library(rgeos)
library(tidyverse)
library(geosphere)

#my points 
points <-read.csv("Analyses/example_points.csv") #this is the table included below of 4 example locations.  

|ID|LON       |LAT       |
|--|----------|----------|
|a |-2.515478 |-69.53887 |
|b |-2.601405 |-69.79783 |
|c |-0.153548 |-69.45126 |
|d |26.06987  |-69.55020 |

#my line
line <- <- readOGR('Environmental_Data/COAST/add_coastline_high_res_polygon_v7_5.shp/') #this is the shapefile linked above

#convert points to spatial object
coordinates(points) <- ~LON+LAT

distance <- geosphere::dist2Line(p = points, line = line, distfun = distGEO)

However, I get an error: "Error in .spDistPoint2Line(p, line, distfun) : Points are projected. They should be in degrees (longitude/latitude)".

The package documentation states that p can be: "longitude/latitude of point(s). Can be a vector of two numbers, a matrix of 2 columns (first one is longitude, second is latitude) or a SpatialPoints object*" - which is what I'm providing it with. I have seen the same issue encountered on a Reddit post (unanswered) but not on here.

My desired output is as below (distances under distance to coast are made up for now!). I have ~3000 locations I need to find the distance to the coastline for.

ID LON LAT Dist_to_coast (km)
a -2.515478 -69.53887 40
b -2.601405 -69.79783 24
c -0.153548 -69.45126 74
d 26.06987 -69.55020 23

Is there an alternative/better means of doing this?

Thank you.


Solution

  • You have loaded sf, any particular reason for not using sf::st_distance() for the task? Would still need to transform though, as there are 4 sample points vs ~140MB shapefile with ~17000 polygons, points were transformed:

    library(ggplot2)
    library(dplyr)
    library(sf)
    
    coastline <- st_read("add_coastline_high_res_polygon_v7_6.shp/")
    p <- readr::read_delim(
     "ID|LON       |LAT      
      a |-2.515478 |-69.53887
      b |-2.601405 |-69.79783
      c |-0.153548 |-69.45126
      d |26.06987  |-69.55020" , delim = "|", trim_ws = T) %>% 
      st_as_sf(coords = c("LON", "LAT"), crs = "WGS84") %>%
      # transform points to match crs of the shapefile
      st_transform(st_crs(coastline))
    
    # number of different surface polygons
    table(coastline$surface)
    #> 
    #>  ice shelf ice tongue       land     rumple 
    #>        325         37      17233         64
    
    # create a single multipolygon, can take a while;
    # you may need to filter first to define any surface types you might want to 
    # include / exclude ("land" also includes islands)
    system.time({
      ucoastline <- st_union(coastline)
    })
    #>    user  system elapsed 
    #>  103.40   11.72  116.08
    
    p$dist_ucoastline <- st_distance(p,ucoastline)
    
    # or perhaps select land polygon with max area to 
    # ignore ice and all the islands:
    land_max <- coastline %>% 
      slice_max(st_area(.))
    p$land_max <- st_distance(p,land_max)
    
    ggplot() +
      geom_sf(data = st_simplify(ucoastline,dTolerance = 1000), fill = "lightblue", color = NA) +
      geom_sf(data = st_simplify(land_max,dTolerance = 1000), fill = "gray70") +
      geom_sf(data = p, shape =4, color="red", size = 5) +
      theme_bw()
    
    

    Result:

    # convert coordinates back to WGS84, 
    # geometries to coordinate columns
    bind_cols(
      st_transform(p, crs = "WGS84") %>% st_coordinates(),
      st_drop_geometry(p)
    ) 
    
    #> # A tibble: 4 × 5
    #>        X     Y ID    dist_ucoastline[,1] land_max[,1]
    #>    <dbl> <dbl> <chr>                 [m]          [m]
    #> 1 -2.52  -69.5 a                  40742.      180479.
    #> 2 -2.60  -69.8 b                  39750.      157043.
    #> 3 -0.154 -69.5 c                   6629.      186878.
    #> 4 26.1   -69.6 d                  45683.      121500.
    
    

    Created on 2022-11-23 with reprex v2.0.2