CODEWITHSUNDEEP
rgeospatialpolygongeosphere
Ellie_Petrels
Ellie_Petrels

Reputation: 17

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.

Upvotes: 0

Views: 331

Answers (1)

margusl
margusl

Reputation: 17389

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

Upvotes: 2

Related Questions