Two methods to get distance from point to enclosing polygon boundary

point_dist_to_polygon.Rmd

---
title: 'C14 analysis: Coastal vs Inland'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Import and clean the data:

```{r import-data}
library(tidyverse)
library(sf)
library(here)

# https://archaeologydataservice.ac.uk/archives/view/austarch_na_2014/
dates <- read_csv(here("analysis", "data", "Austarch_1-3_and_IDASQ_28Nov13-1.csv"))

# fix ages,  remove empty coords, convert to spatial
dates_sf <- 
  dates %>% 
  mutate(AGE = parse_number(AGE)) %>% 
  filter(!is.na(LONGITUDE)) %>% 
  st_as_sf(coords = c('LONGITUDE', 'LATITUDE'),
                     remove = FALSE)
  

# from https://www.naturalearthdata.com/downloads/10m-physical-vectors/
coastline <- st_read(here("analysis", "data", "ne_10m_coastline/ne_10m_coastline.shp"))

# trim world shapefile so it's only Australia
lower_left <-  c(-45, 110)
upper_right <- c(-10, 159)
oz_coastline <- 
  st_crop(coastline, 
          c(xmin = lower_left[2], xmax = upper_right[2], 
            ymin = lower_left[1], ymax = upper_right[1]))

# check it
ggplot() +
  geom_sf(data = oz_coastline) 

# simplify and exclude islands
oz_coastline <- 
  oz_coastline %>% 
  slice(11:14) %>% 
  st_combine %>% 
  st_union 

# simplify further
oz_coastline <- 
  rmapshaper::ms_simplify(as(oz_coastline, "Spatial"), 
                          keep_shapes = TRUE, 
                          keep = 0.01) %>% 
  st_as_sfc()

# check it
ggplot() +
  geom_sf(data = oz_coastline) 

# we need to set the Coordinate Reference System to be the same
st_crs(dates_sf) <- st_crs(oz_coastline)

```

Draw a map:

```{r}
# map coast
ggplot() +
  geom_sf(data = oz_coastline) +
  geom_sf(data = dates_sf,
          aes(colour = AGE_NORM)) +
  coord_sf() +
  theme_minimal() +
  scale_color_viridis_c()
  
```


Get distance of sites from coast:



```{r}

#--------------------------------------

points_sf <- sample_n(dates_sf, nrow(dates_sf))

# get only the lines that are the coast (no islands),
# and combine into one multilinestring
polygon_sf <- 
  oz_coastline 

# very fast!
dist <- st_distance(points_sf, 
                    polygon_sf)

# convert the distance to numeric
points_sf_dist <- 
points_sf %>% 
  mutate(dist = as.numeric(dist) / 1000)

# get lines to show distances
dist_lines <- st_nearest_points(points_sf, polygon_sf)

# visual check that distance-to-coast is 
# computed as expected
ggplot() +
 geom_sf(data = points_sf_dist,
         aes(colour = dist)) +
 geom_sf(data = polygon_sf)  +
 geom_sf(data = dist_lines,
         aes(colour = points_sf_dist$dist)) +
 scale_color_viridis_c() +
 coord_sf() +
 theme_minimal() 


#--------------------------------------

gdist <- 
  geosphere::dist2Line(p = st_coordinates(points_sf), 
                           line =  polygon_sf %>% 
                                    st_coordinates %>% 
                                    as_tibble %>% 
                                    select(X, Y))


# combine initial site point data with distance to coastline
df <- bind_cols(points_sf, as.data.frame(gdist))

# take a look
ggplot() + 
  geom_sf(data = polygon_sf) +
  geom_sf(data = points_sf) +
  geom_segment(data = df,
               aes(x = LONGITUDE,
                   y = LATITUDE,
                   xend = lon,
                   yend = lat,
                   colour = distance))



#--------------------------------------


```