Antimeridian problem

CHANGELOG.md

@@ -7,6 +7,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/openghg/openghg/compare/0.1.1...HEAD)
 
+### Fixed
+- Added a function to ensure correct polygon display across the dateline - [PR #285](https://github.com/4DModeller/fdmr/pull/285)
+
 ### Added
 
 - Added mouse pointer coordinates header and standard measurement tool - [PR #260](https://github.com/4DModeller/fdmr/pull/260)

DESCRIPTION

@@ -58,7 +58,8 @@ Imports:
     bslib,
     bsicons,
     future,
-    leafem
+    leafem,
+    sfheaders
 Suggests: 
     bookdown,
     knitr,

NAMESPACE

@@ -1,6 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
+export(antimeridian_wrapping)
 export(clear_caches)
+export(convert_from_lon_360)
 export(create_prediction_field)
 export(get_tutorial_datapath)
 export(latlong_to_utm)

R/plot_mapping.R

@@ -12,6 +12,7 @@
 #' @param polygon_line_colour Polygon surrounding line colour
 #' @param polygon_line_weight Polygon surrounding line weight
 #' @param reverse Reverse the colour palette if TRUE
+#' @param wrapping Split polygons along the antimeridian (-180/180 boundary) if TRUE
 #'
 #' @return leaflet::leaflet
 #' @export
@@ -26,7 +27,8 @@ plot_map <- function(polygon_data = NULL,
                      polygon_line_colour = "grey",
                      polygon_line_weight = 1,
                      polygon_fill_opacity = 0.6,
-                     reverse = FALSE) {
+                     reverse = FALSE,
+                     wrapping = FALSE) {
   if (is.null(polygon_data) && is.null(raster_data)) {
     stop("Polygon or raster data must be given.")
   }
@@ -41,6 +43,9 @@ plot_map <- function(polygon_data = NULL,
   layers <- c()
 
   if (!is.null(polygon_data)) {
+    if (isTRUE(wrapping)) {
+      polygon_data <- fdmr::antimeridian_wrapping(polygon_data, crs = "+proj=longlat +datum=WGS84", unique_inst = TRUE, to_sp = FALSE)
+      }
     if (!is.null(domain)) {
       colours <- leaflet::colorNumeric(palette = palette, domain = domain, reverse = reverse)
       polygon_fill_colour <- ~ colours(domain)

R/plot_mesh.R

@@ -24,7 +24,9 @@ plot_mesh <- function(mesh, spatial_data = NULL, longitude_column = "LONG", lati
     )
   }
 
-  spatial_mesh <- fdmr::mesh_to_spatial(mesh = mesh, crs = expected_crs)
+  spatial_mesh_original <- fdmr::mesh_to_spatial(mesh = mesh, crs = expected_crs)
+
+  spatial_mesh <- fdmr::antimeridian_wrapping(spatial_mesh_original, crs = expected_crs, unique_inst = FALSE, to_sp = FALSE)
 
   plot_polygons <- FALSE
   plot_points <- FALSE

R/shiny_meshbuilder.R

@@ -174,7 +174,7 @@ meshbuilder_shiny <- function(
     mesh_spatial <- shiny::reactive(
       suppressMessages(
         suppressWarnings(
-          fdmr::mesh_to_spatial(mesh = mesh(), crs = crs)
+          fdmr::antimeridian_wrapping(fdmr::mesh_to_spatial(mesh = mesh(), crs = crs), crs = crs, unique_inst = FALSE, to_sp = FALSE)
         )
       )
     )
@@ -183,7 +183,7 @@ meshbuilder_shiny <- function(
       m <- leaflet::leaflet()
       m <- leaflet::addTiles(m, group = "OSM")
       m <- leaflet::addPolygons(m, data = mesh_spatial(), weight = 0.5, fillOpacity = 0.2, fillColor = "#5252ca", group = "Mesh")
-      m <- leaflet::addMeasure(position = 'bottomleft', primaryLengthUnit = 'kilometers', primaryAreaUnit = 'sqmeters')
+      m <- leaflet::addMeasure(m, position = 'bottomleft', primaryLengthUnit = 'kilometers', primaryAreaUnit = 'sqmeters')
       m <- leafem::addMouseCoordinates(m, native.crs = TRUE)
       if (plot_polygons) {
         m <- leaflet::addPolygons(m, data = spatial_data, fillColor = "#d66363", color = "green", weight = 1, group = "Spatial")

R/util_coords.R

@@ -40,3 +40,112 @@ has_coords <- function(spatial_data) {
     }
   )
 }
+
+#' Convert longitudes from 0 to 360 degrees to -180 to 180 degrees
+#'
+#' @param sf_data An sf object; does not accept SpatialPolygon* objects
+#' @param crs CRS as a proj4string or EPSG code
+#' @param add_data Select if data associated with the object are carried forward by the transformed version, defaults to FALSE
+#' @param longitude_column Name of longitude, defaults to LONG
+#'
+#' @return polygons with converted coordinates
+#' @export
+convert_from_lon_360 <- function(sf_data, crs = 4326, add_data = TRUE, longitude_column = "LONG") {
+  # Get polygon coordinates
+  coords <- sf::st_coordinates(sf_data)
+
+  # Check if POLYGON contains holes, remove if so
+  if ((base::length(base::unique(coords[, "L1"])) > 1) & (base::ncol(coords) > 3)) {
+    warning("Converter does not define the interior holes. Please convert your data yourself if these features are crucial")
+    sf_data <- sfheaders::sf_remove_holes(obj = sf_data)
+    coords <- sf::st_coordinates(sf_data)
+  }
+
+  coords_conv <- coords
+
+  # Convert data coordinates if selected
+  if (base::isTRUE(add_data) & (!base::is.null(longitude_column))) {
+    if (base::is.null(sf_data[["LONG"]]) & (longitude_column == "LONG")) {
+      warning("`LONG` not found. Please change the default `longitude_column` name to the longitude name in your dataset")
+    } else {
+      sf_data[[longitude_column]] <- base::ifelse(sf_data[[longitude_column]] > 180, 
+                                                  sf_data[[longitude_column]] - 360, 
+                                                  sf_data[[longitude_column]])
+    }
+  }
+
+  # Convert polygon longitudes
+  coords_conv[, "X"] <- base::ifelse(coords[, "X"] > 180, coords[, "X"] - 360, coords[, "X"])
+
+  # Define the feature column and convert to POLYGON
+  if (base::ncol(coords_conv) == 4) {
+    poly <- "L2"
+    conv <- sfheaders::sf_polygon(obj = coords_conv, x = 'X', y = 'Y', polygon_id = poly)
+    conv <- sf::st_sf(conv, crs = crs)
+  } else if (base::ncol(coords_conv) == 5) {
+    poly <- "L3"
+    multipoly <- "L2"
+    conv <- sfheaders::sf_multipolygon(obj = coords_conv, x = 'X', y = 'Y', multipolygon_id = multipoly, polygon_id = poly)
+    conv <- sf::st_sf(conv, crs = crs)
+  } else {
+    stop("Object type not supported (must be POLYGON or MULTIPOLYGON)")
+  }
+
+  # Add the associated data if requested
+  if (base::isTRUE(add_data)) {
+    r_conv <- sf_data
+    r_conv$geometry <- conv$geometry
+  } else {
+    r_conv <- conv
+  }
+
+  # Return a list of POLYGON objects with shifted longitudes
+  return(r_conv)
+}
+
+#' Split polygons into two if crossing the dateline (antimeridian)
+#' 
+#' @param polygon_data SpatialPolygon or SpatialPolygonDataFrame object
+#' @param crs CRS as a proj4string or EPSG code, default EPSG:4326 (datum=WGS84)
+#' @param unique_inst Option to remove duplicate polygons (with the same coordinates), default TRUE
+#' @param to_sp Option to convert POLYGON object into SpatialPolygon*, default TRUE
+
+#' @return SpatialPolygon or SpatialPolygonDataFrame split along the antimeridian
+#' @export
+antimeridian_wrapping <- function(polygon_data, crs = 4326, unique_inst = TRUE, to_sp = TRUE) {
+  # Convert to sf object if needed
+  if (!is(polygon_data, 'sf')) {
+    if (base::isTRUE(unique_inst)) {
+      all_coords <- base::lapply(polygon_data@polygons, function(x) {
+        base::lapply(x@Polygons, function(x) {
+          sp::coordinates(x)
+        })
+      })
+      polygon_data <- polygon_data[!base::duplicated(all_coords),]
+    }
+    sf_data <- sf::st_as_sf(polygon_data)
+    if (!isTRUE(sf::st_is_longlat(sf_data))) {
+      stop("Projected coordinates not allowed, please choose a different projection to avoid the antimeridian problem")
+    }
+  } else {
+    sf_data <- polygon_data
+  }
+
+  # Check if longitudes go from -180 to 180 rather than 0 to 360
+  if (sf::st_bbox(sf_data)[['xmin']] > 0) {
+    polygon_data_conv <- convert_from_lon_360(sf_data, crs = crs, add_data = FALSE)
+    warning("Polygon coordinates [0;360] have been converted to [-180;180]")
+  }
+
+  # Apply the function
+  sf_split <- sf::st_wrap_dateline(sf_data)
+
+  # Revert sf objects back to sp
+  if (base::isTRUE(to_sp)) {
+    polygon_data_conv <- sf::as_Spatial(sf_split)
+  } else {
+    polygon_data_conv <- sf_split
+  }
+
+  return(polygon_data_conv)
+}