calculate distance of points to polygon boundary using terra package in R-CodePudding

I am trying to calculate distance of points within a country to country boundary

library(terra)
library(geodata)
library(ggplot2)
library(geodata)

# get a shapefile of a country 
gabon <- geodata::gadm('GAB', level = 0, path = getwd())    
canvas <- terra::rast(xmin = ext(gabon)[1], 
                      xmax = ext(gabon)[2], 
                      ymin = ext(gabon)[3], 
                      ymax = ext(gabon)[4],
                      resolution = 0.08,
                      crs = crs(gabon),
                      vals = 0)
pts <- as.points(canvas)    
pts <- terra::crop(pts, gabon) # extract the points in the limits of Gabon    

plot(pts)
plot(gabon, border = "blue", add = T)

I want to calculate shortest distance of each point in pts to the boundary of the country

gabon_lines <- terra::as.lines(gabon)

# calculation of the distance between the boundary and points
dis_pts <- terra::distance(pts, gabon_lines, pairwise = FALSE, unit="km")
range(dis_pts)
# 0.00000046 1.63706213. seems quite low 

dat <- data.frame(dist = as.vector(dis_pts), 
                  crds(pts))

col_dist <- brewer.pal(11, "RdGy")

ggplot(dat, aes(x, y, fill = dist))   #variables
  geom_tile()   #geometry
  scale_fill_gradientn(colours = rev(col_dist))  # colors for plotting the distance
  labs(fill = "Distance (km)")  #legend name
  theme_void()  #map theme
  theme(legend.position = "bottom") #legend position

I think the range of distance I am getting is very low since Gabon is quite big so I was expecting distance of points in the middle to be larger. Is there anything I am doing wrong here?

CodePudding user response：

The problem seems to be with the crs used. The result you have above is accurate, but the units are in degrees (latitude & longitude). A relatively quick fix is to reproject the data using

The dimensions come out a little wonky since the plot isn't using a crs. Changing the data to sf points makes things look a little better:

library(sf)

st_as_sf(dat, coords = c("x", "y")) %>% 
  ggplot()   
    geom_sf(aes(color = dist))   
    scale_color_gradientn(colours = rev(col_dist))

CodePudding user response：

That needs to be fixed, but you can do this

library(terra)
library(geodata)

# get a shapefile of a country 
gabon <- geodata::gadm('GAB', level = 0, path = getwd())    
canvas <- terra::rast(gabon, resolution = 0.08, vals=0)
m <- mask(canvas, gabon, inverse=TRUE)
d <- distance(m)
plot(d)