CODEWITHSUNDEEP
rr-sf
James S
James S

Reputation: 25

R: st_buffer is giving a wildly inaccurate result when projected to CRS WSG84

I have a latitude and longitude that I am treating as a center point. I want to draw a 100 nautical mile circle around that point. Since 1 degree of lat/lon is approximately 60 nautical miles, I would expect the difference between the max and min latitude (as well as the max and min longitudes) to be a little over 3 degrees.

When I create a 185,200 meter (1,852 m per nautical mile) buffer around the lat/lon (28.5, 86.5), the bounding box of the resulting polygon has a min/max lon of 84.85 and 88.16, respectively. This makes sense. However, the resulting min/max lat is 0.24 and 56.78. That makes no sense at all to me and I have no idea what's going on.

This is the code I've tried

test_point = st_as_sf(data.frame(Lat=28.5,Lon=86.5),coords=c("Lat","Lon"),remove=FALSE)
st_crs(test_point) = 4326
buff1 = st_buffer(test_point,100*1852)
buff2 = st_transform(buff1,crs=4326)

When I look at the min/max x,y coordinates for buff1, the differences in the min/max between each dimension is about 200 nautical miles; this makes sense since that's twice the radius I've specified. But things go haywire when I go back to WSG84

buff1
#Simple feature collection with 1 feature and 2 fields
#Geometry type: POLYGON
#Dimension:     XY
#Bounding box:  xmin: 667639.2 ymin: 9981165 xmax: 1038039 ymax: 10351570
#CRS:           +proj=utm +zone=16
#   Lat  Lon                       geometry
#1 28.5 86.5 POLYGON ((1038039 10166365,...

#(Note that (10351570-9981165)/1852 is about 200 and (1038039-667639.2)/1852 is about 200 too.

buff2
#Simple feature collection with 1 feature and 2 fields
#Geometry type: POLYGON
#Dimension:     XY
#Bounding box:  xmin: 0.2407929 ymin: 84.84537 xmax: 56.77518 ymax: 88.15623
#Geodetic CRS:  WGS 84
#   Lat  Lon                       geometry
#1 28.5 86.5 POLYGON ((20.35843 84.95627...

# Note how far apart the xmin and xmax are.

Upvotes: 0

Views: 219

Answers (1)

margusl
margusl

Reputation: 17399

That test_point ends up near North Pole, at 86N, and up there circle with 100 nmile radius can indeed cover 57+ degrees longitude (black points are bbox corners):

library(sf)
library(ggplot2)
library(mapview)

test_point = st_as_sf(data.frame(Lat=28.5,Lon=86.5),coords=c("Lat","Lon"),remove=FALSE)
st_crs(test_point) = 4326
buff1 = st_buffer(test_point,100*1852)
buff1
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 0 ymin: 84.81199 xmax: 57.29792 ymax: 88.18218
#> Geodetic CRS:  WGS 84
#>    Lat  Lon                       geometry
#> 1 28.5 86.5 POLYGON ((2.271324 86.23322...

# let's visualize
bbox_points <- 
  st_bbox(buff1) %>% 
  st_as_sfc() %>% 
  st_cast("POINT") %>% 
  st_as_sf() %>% 
  mutate(coord = lapply(x, \(x) paste0(round(x[1],3),"\n", round(x[2],3))))

ggplot() +
  geom_sf(data = buff1, alpha = .2) +
  geom_sf(data = test_point, color = "red") +
  geom_sf(data = bbox_points) +
  geom_sf_text( data = bbox_points, aes(label = coord)) +
  coord_sf(crs = 6115) +
  theme_bw() +
  theme(panel.border = element_blank())

Assuming from UTM Zone16 CRS in your code output, you might have meant 86.5 West not 86.5 East for longitude, so after fixing coordinate order (coords parameter in st_as_sf() uses x,y order) and changing sign for longitude, we now end up with this:

test_point_2 <- st_as_sf(
  data.frame(Lat =  28.5,
             Lon = -86.5),
  coords = c("Lon","Lat"),
  remove = FALSE,
  crs = "WGS84")

buff_2 <- st_buffer(test_point_2, units::set_units(100, nmile))
buff_2
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -88.4016 ymin: 26.81976 xmax: -84.5851 ymax: 30.18484
#> Geodetic CRS:  WGS 84
#>    Lat   Lon                       geometry
#> 1 28.5 -86.5 POLYGON ((-84.94013 29.4612...

# visualize
buff_extend <- st_buffer(test_point_2, units::set_units(500, nmile))
viewExtent(buff_extend) +
mapview(buff_2) +
  mapview(test_point_2)

x-range of the bounding box is much closer to your expected result now, though 1 degree of lat being ~60 nautical miles only works at the equator.

Created on 2023-03-11 with reprex v2.0.2

Upvotes: 1

Related Questions