CODEWITHSUNDEEP
rmapsgeospatialr-sf
max946
max946

Reputation: 33

Projecting a quartered circle with a 50km radius in r/sf?

I'm hoping to create a series of quartered circles (i.e. circles split into 4 equal quadrants), each with a 50km radius, that I can map onto various longitudes and latitudes throughout the United States. I'd also like the option to rotate these quartered circles as desired.

Using the code below (and guidance from here), I've been able to make the following start:

New York State Map

I have two questions:

  1. How can I meaningfully set the radius of these circles? Is there a way to draw shapes a certain distance (in km) from a coordinate in a projected CRS? So far I'm defining the radius in terms of degrees of longitude and latitude, but distance would be more useful.
  2. My circles appear to be turning into ellipses after projecting them and mapping them in WGS84. Is there any way to prevent this from happening?

I would be happy to consider alternative approaches. Thanks!

library(sf)
library(ggplot2)
library(maps)

#Two functions to create coordinate quartered circle polygons

#x = long, y = lay, r = radius, theta_rotate = rotation
st_wedge <- function(x,y,r,start,width, theta_rotate){
   n <- 20
      theta = seq(start+theta_rotate, start+width+theta_rotate, length=n)
      xarc = x + r*sin(theta) 
      yarc = y + r*cos(theta)
      xc = c(x, xarc, x) 
      yc = c(y, yarc, y)
      st_polygon(list(cbind(xc,yc)))   
}


st_wedges <- function(x, y, r, nsegs, theta_rotatex){
   width = (2*pi)/nsegs
   starts = (1:nsegs)*width
   polys = lapply(starts, function(s){st_wedge(x,y,r,s,width, theta_rotatex)})
   #Cast to crs 4326, WGS84
   mpoly = st_cast((st_sfc(polys, crs = 4326)), "MULTIPOLYGON")
   mpoly
}

#Create quartered sf circle polygon 
custom_circle_sf <- st_wedges(x = -76, y = 43, r = .3, nsegs = 4, theta_rotatex = 200) %>%
 st_sf() %>% 
   mutate(group = row_number()) %>% dplyr::select(group, geometry)

#Create New York State sf polygon
ny_map_sf <- map_data("state", region="new york")  %>% 
st_as_sf(coords = c("long", "lat"), crs = 4326)   %>% 
   group_by(group) %>%
   summarise(geometry = st_combine(geometry)) %>% 
   st_cast("POLYGON")

#Plot results
ggplot() +
   geom_sf(data=ny_map_sf,
           size = 1,
           colour = "blue", 
           fill   = "white") + 
   geom_sf(data=custom_circle_sf,
           size = .1,
           aes(fill=group),
           colour = "white")

Upvotes: 1

Views: 347

Answers (2)

max946
max946

Reputation: 33

For anyone who is curious about splitting polygons in sf using R, this was how I went about solving this:

#Function to create circle with quadrants. Save desired projection as projected_crs
create_circle <- function(lat_x, long_y, theta_x, buffer_m){
   #Create circle with radius buffer_m centered at (lat_x, long_y)
   circle_buffer  <- st_point(c(lat_x, long_y)) %>%  st_sfc(crs = 4326) %>% 
      st_cast("POINT")  %>% 
      st_transform(projected_crs) %>%
      st_buffer(buffer_m)
   
   #Create two orthogonal lines at origin 
   p1 <- rbind(c(lat_x,long_y - 1), c(lat_x,long_y + 1))
   p2 <- rbind(c(lat_x+1,long_y), c(lat_x-1,long_y))
   mls <- st_multilinestring(list(p1,p2))  %>% st_sfc(crs = 4326) %>% 
      st_transform(projected_crs) 
   
   #Use orthogonal lines to split circle into 4 quadrants
   x1 <- st_split(circle_buffer, mls) 
   
   #Convert origin into projected CRS
   center_in_crs  <- st_point(c(lat_x, long_y)) %>% 
      st_sfc(crs = 4326) %>%
      st_transform(projected_crs)
   
   sp_obj <- x1 %>% st_collection_extract(type="POLYGON") %>%
      #Convert to spatial to use sp functions
      as_Spatial() %>% 
      #rotate x degrees
      elide(rotate = theta_x + 45, center  = center_in_crs[[1]]) %>% 
      #return to sf 
      st_as_sf()

Upvotes: 1

Micha
Micha

Reputation: 493

Regarding your question 2: "circles appear to be turning into ellipses". If you add to your ggplot the coord_equal() function then the grid will be square, and the ellipses will be shown as circles.

Upvotes: 0

Related Questions