CODEWITHSUNDEEP
rgisr-sf
Nicolas Roelandt
Nicolas Roelandt

Reputation: 151

Rstats - How to convert kilometers into arc degrees in order to create buffers with {sf} / {units}?

I have a csv file with cities and their coordinates in WGS84. I would like to create buffers of 20 kilometers around in order to find too close cities. In the example below, Islamabad and Rawalpindi are too close to each other.

I was able to create the df and the geom but when I call st_buffer(), it tells me that it cannot convert km into degrees angle.

I tried using the units package but it doesn't seem to handle degrees angles.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(units)
#> udunits system database from D:/Roelandt/Documents/R/win-library/3.5/units/share/udunits

# Create a tribble from the data

df <- tibble::tribble(
   ~ID,             ~city,     ~lat,    ~lon,
  1172,           "Zaria", 11.11128,   7.7227,
  1173,            "Oslo", 59.91273, 10.74609,
  1174, "Masqat (Muscat)", 23.61387,  58.5922,
  1175,      "Bahawalpur",     29.4, 71.68333,
  1181,"Islamabad",33.70351,73.059373,
  1194,"Rawalpindi",33.6,73.0666667

  )
df
#> # A tibble: 6 x 4
#>      ID city              lat   lon
#>   <dbl> <chr>           <dbl> <dbl>
#> 1  1172 Zaria            11.1  7.72
#> 2  1173 Oslo             59.9 10.7 
#> 3  1174 Masqat (Muscat)  23.6 58.6 
#> 4  1175 Bahawalpur       29.4 71.7 
#> 5  1181 Islamabad        33.7 73.1 
#> 6  1194 Rawalpindi       33.6 73.1

cities_df = st_as_sf(df, coords = c("lon", "lat"), crs = 4326)

# buffer size 
a = as_units(20, "km")

#Create buffers
cities_buffers <- cities_df %>% 
  st_buffer(dist = a)
#> Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
#> endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
#> Error: cannot convert km into °
cities_buffers
#> Error in eval(expr, envir, enclos): objet 'cities_buffers' introuvable

Created on 2019-01-21 by the reprex package (v0.2.1)

I would like to know if there is a way to convert the buffer size in geographical degrees.

I was thinking about converting the dataset to a projected crs like UTM, but I'm not sure it is a good idea on worldwide data.

Thanks,

Nicolas

Upvotes: 4

Views: 1133

Answers (1)

sebdalgarno
sebdalgarno

Reputation: 3197

From Jessie Sadler's great blog post (https://www.jessesadler.com/post/simple-feature-objects/): "Beginning with version 0.6 of the sf package, st_distance() uses the lwgeom package, which in turn links to geometric functions from the liblwgeom library used by PostGIS, to make geometric calculations on longitude and latitude values."

I would approach this by getting a distance matrix rather than creating buffers.

dist_mat <- st_distance(cities_df)
#> dist_mat
#Units: [m]
#        [,1]    [,2]    [,3]      [,4]       [,5]       [,6]
#[1,]       0 5421407 5550432 6898123.8 7057404.43 7057248.70
#[2,] 5421407       0 5461507 5616666.0 5306301.50 5315880.95
#[3,] 5550432 5461507       0 1452564.2 1799547.86 1793603.09
#[4,] 6898124 5616666 1452564       0.0  494716.65  483860.32
#[5,] 7057404 5306302 1799548  494716.7       0.00   11500.84
#[6,] 7057249 5315881 1793603  483860.3   11500.84       0.00

From the matrix, you can then determine which points are <20,000m apart

> dist_mat < set_units(20000, "m") & dist_mat > set_units(0, "m")
#      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
#[1,] FALSE FALSE FALSE FALSE FALSE FALSE
#[2,] FALSE FALSE FALSE FALSE FALSE FALSE
#[3,] FALSE FALSE FALSE FALSE FALSE FALSE
#[4,] FALSE FALSE FALSE FALSE FALSE FALSE
#[5,] FALSE FALSE FALSE FALSE FALSE  TRUE
#[6,] FALSE FALSE FALSE FALSE  TRUE FALSE

Upvotes: 2

Related Questions