Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve check_antimeridian #28

Open
jflowernet opened this issue Dec 19, 2024 · 0 comments
Open

Improve check_antimeridian #28

jflowernet opened this issue Dec 19, 2024 · 0 comments

Comments

@jflowernet
Copy link
Member

#check_antimeridian incorrectly determines that bounding box crossing the antimeridian does not

kiribati <- spatialgridr::get_boundary(name = "Kiribati", country_type = "sovereign")

kir_crs <- '+proj=laea +lon_0=-159.609375 +lat_0=0 +datum=WGS84 +units=m +no_defs'

kir_local_crs <- sf::st_transform(kiribati, kir_crs)

kir_local_bbox <- sf::st_bbox(kir_local_crs) |>
  sf::st_bbox() |>
  sf::st_as_sfc() |>
  sf::st_as_sf()

#returns FALSE, i.e. does not cross antimeridian (but it does!)
spatialgridr:::check_antimeridian(kir_local_bbox, 4326)
#> [1] FALSE

#check_antimeridian projects the bounding box of the polygon into epsg:4326, but results in a bounding box that is 'inverse' - this is the code it uses:
sf::sf_use_s2(TRUE)
sf::st_transform(kir_local_bbox, 4326) |> 
  sf::st_bbox() |>
  sf::st_as_sfc() |>
  plot(axes = T)

sf::st_transform(kir_local_bbox, 4326) |> 
  sf::st_bbox()
#>        xmin        ymin        xmax        ymax 
#> -146.868970  -13.789183  167.699940    7.830588

#this behaviour is only for the bbox, the polygon boundaries work ok

spatialgridr:::check_antimeridian(kir_local_crs, 4326)
#> [1] TRUE

sf::st_transform(kir_local_crs, 4326) |>
  sf::st_bbox()
#>        xmin        ymin        xmax        ymax 
#> -180.000000  -13.838333  179.998831    7.879056

#the problem is likely due to lack of nodes - when polygon is split at antimeridian, it is broken into two sets of points either side of the antimeridian which are then re-joined after projecting, but the join is the 'wrong' way since it cannot be joined across the antimeridian in 4326

#try densifying before transforming - this add nodes and fixes the problem
kir_local_bbox_dense <- terra::densify(terra::vect(kir_local_bbox), interval = 1e4) |>
  sf::st_as_sf()

spatialgridr:::check_antimeridian(kir_local_bbox_dense, 4326)
#> [1] TRUE

sf::st_transform(kir_local_bbox_dense, 4326) |> 
  sf::st_bbox()
#>        xmin        ymin        xmax        ymax 
#> -179.962616  -13.876564  179.969959    7.879078

sf::st_transform(kir_local_bbox_dense, 4326) |> 
  plot(axes = T)

#could improve check_antimeridian code using this idea: https://gis.stackexchange.com/questions/256329/determining-whether-polygon-crosses-the-antimeridian

Created on 2024-12-19 with reprex v2.1.1

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant