0

I'm trying to create a polygon that has longitudinal limits as 150, -170, i.e. crosses the 180 meridian dateline.

I've tried:

x = c(-170, -170, 150, 150) #long limits

y = c(-25,-57,-57,-25) #lat limits

polygon = cbind(x, y) %>%
  st_linestring() %>%
  st_cast("POLYGON") %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES")) %>%   #thought this line could solve it
  st_sfc(crs = 4326, check_ring_dir = TRUE) %>%
  st_sf()

That is not solving the problem, even if I delete the 'st_cast' cast line or use 'MULTIPOLYGONS' instead of 'POLYGONS' in there. I've also created one polygon with positive longitudes and another for the negative ones, and then combined them, but that's not working well (R runs it, but I get nothing when plotting the object with combined polygons).

I would greatly appreciate it if you could provide your ideas on this :)

dieghernan
  • 2,690
  • 8
  • 16
Elisa
  • 1
  • 1
  • Ok...it worked when I used 'st_union' to combine the polygons with positive and negative longitudes (I was trying 'st_intersection' before). But if someone knows how to do it without the need to create and combine polygons, that would be great to learn. – Elisa Mar 28 '22 at 01:40

1 Answers1

2

I think you can do it simply passing the coordinates, but this should be combined also with the right coordinate reference system:

  • You can create a POLYGON from a bounding box quickly, as far as you assign it the class of a bounding box box and after that use st_as_sfc().
  • Use sf_use_s2(TRUE), available on sf >= 1.0.0.

See here how to do it:

library(sf)
sf_use_s2(TRUE)

# From bounding box
box <- c(xmin=-170, ymin=-57, xmax=150, ymax=-25)
class(box) <- "bbox"

box_end <- box |>  
  st_as_sfc() |> 
  st_as_sf(crs=4326) 

box_end
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -170 ymin: -57 xmax: 150 ymax: -25
#> Geodetic CRS:  WGS 84
#>                                x
#> 1 POLYGON ((-170 -57, 150 -57...

# Check if point in polygon

# This point should not be on your polygon
ptest1 <- st_as_sfc("POINT(130 -40)", crs=4326)
st_contains(box_end, ptest1) 
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `contains'
#>  1: (empty)

# This should be on the polygon
ptest2 <- st_as_sfc("POINT(-176 -40)", crs=4326)
st_contains(box_end, ptest2) 
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `contains'
#>  1: 1

If you want to plot it you must use a suitable projection for your coordinates. In this case I use an Ortographic projection centered in c(-160, -40):

# Just for example: Using Pacific centered crs
library(ggplot2)
library(giscoR)

data("gisco_countries")

ggplot(gisco_countries) +
  geom_sf()+
  geom_sf(data=box_end, fill="red") +
  geom_sf(data=ptest1, col="green", size=2) +
  geom_sf(data=ptest2, col="blue", size=2) +
  coord_sf(crs = "+proj=ortho +x_0=0 +y_0=0 +lat_0=-40 +lon_0=160")

enter image description here

If you want to have a MULTIPOLYGON instead a POLYGON

Use st_shift_longitude() + st_wrap_dateline().

library(sf)
sf_use_s2(TRUE)

# From bounding box
box <- c(xmin=-170, ymin=-57, xmax=150, ymax=-25)
class(box) <- "bbox"

box_end <- box |>  
  st_as_sfc() |> 
  st_as_sf(crs=4326) |> 
  # This splits the POLYGON and creates a MULTIPOLYGON
  # Note that the bounding box is also affected
  st_shift_longitude() |> 
  st_wrap_dateline()

box_end
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -180 ymin: -57 xmax: 180 ymax: -25
#> Geodetic CRS:  WGS 84
#>                               x
#> 1 MULTIPOLYGON (((150 -57, 15...

# Check if point in polygon

# This point should not be on your polygon
ptest1 <- st_as_sfc("POINT(130 -40)", crs=4326)
st_contains(box_end, ptest1) 
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `contains'
#>  1: (empty)

# This should be on the polygon
ptest2 <- st_as_sfc("POINT(-176 -40)", crs=4326)
st_contains(box_end, ptest2) 
#> Sparse geometry binary predicate list of length 1, where the predicate
#> was `contains'
#>  1: 1

# Just for example: Using Robinson
library(ggplot2)
library(giscoR)

data("gisco_countries")

ggplot(gisco_countries) +
  geom_sf()+
  geom_sf(data=box_end, fill="red") +
  geom_sf(data=ptest1, col="green", size=2) +
  geom_sf(data=ptest2, col="blue", size=2) +
  coord_sf(crs = "+proj=robin")

enter image description here

dieghernan
  • 2,690
  • 8
  • 16