CODEWITHSUNDEEP
rrasterr-sfterra
89_Simple
89_Simple

Reputation: 3805

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

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)

enter image description here

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

enter image description here

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?

Upvotes: 1

Views: 1248

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47071

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)

Upvotes: 4

mrhellmann
mrhellmann

Reputation: 5489

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 crs 5223.

Most of the code below is copied, changes are below ####

library(terra)
library(geodata)
library(ggplot2)
library(scales)
library(RColorBrewer)

# 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)    

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


#### 
# reproject pts & gabon lines to this new crs:
new_crs <- "+proj=tmerc +lat_0=0 +lon_0=12 +k=0.9996 +x_0=500000 +y_0=500000 +datum=WGS84 +units=m +no_defs +type=crs"

pts2 <- terra::project(pts, new_crs)
gabon_lines2 <- terra::project(gabon_lines, new_crs)

# calculation of the distance between the boundary and points
dis_pts <- terra::distance(pts2, gabon_lines2, pairwise = FALSE, unit="km")
range(dis_pts)
## Now from 1 to about 180 km
## a quick check on google maps & the interior of Gabon is ~180km from the nearest border

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

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

## Not much change from the plot before, but lat & lon degrees are approximately the same near the equator
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

enter image description here

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))

enter image description here

Upvotes: 1

Related Questions