I am using the R package sf. I have a bounding box for Australia's EEZ in Mollweide projection, which I make into an sf polygon and then re-project to long-lat (EPSG 4326) coordinate system. This results in a triangle polygon, as shown in the following reprex. I understand that re-projecting will distort the polygon from its rectangular shape, but I don't understand why it results in a triangle?
Looking at the 2 polygons' coordinates using st_coordinates
I can see that the polygon in lon-lat has only 4 coordinates, while that in Mollweide projection has 5; but I don't understand why one coordinate has been dropped.
Any help appreciated.
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE
library(rnaturalearth)
#> Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
#get australia country polygon in lon lat
aus_lonlat <- ne_countries(country = "australia", returnclass = "sf") |>
st_geometry()
#project to Mollweide
aus_moll <- aus_lonlat |>
st_transform("ESRI:54009")
#define bounding box I'm interested in - this is for Australia's EEZ and is in Mollweide projection
polygon_bbox_moll <- st_bbox(c(xmin =9797297, xmax = 15540552, ymin = -5576498, ymax = -1096691), crs = st_crs("ESRI:54009")) |>
st_as_sfc() |>
st_as_sf()
#plot everything - Mollweide
plot(polygon_bbox_moll, lty = "dotted", axes=T)
plot(aus_moll, add=T)
#project bounding box to lon lat coordinates
polygon_bbox_lonlat <- polygon_bbox_moll|>
st_transform(4326)
#plot everything - lon-lat
plot(polygon_bbox_lonlat, lty = "dotted", axes=T)
plot(aus_lonlat, add=T)
#lon lat polygon has one less coordinate
st_coordinates(polygon_bbox_moll)
#> X Y L1 L2
#> [1,] 9797297 -5576498 1 1
#> [2,] 15540552 -5576498 1 1
#> [3,] 15540552 -1096691 1 1
#> [4,] 9797297 -1096691 1 1
#> [5,] 9797297 -5576498 1 1
st_coordinates(polygon_bbox_lonlat)
#> X Y L1 L2
#> [1,] 124.37141 -47.193643 1 1
#> [2,] 156.21908 -8.883331 1 1
#> [3,] 98.48587 -8.883331 1 1
#> [4,] 124.37141 -47.193643 1 1
<sup>Created on 2023-08-24 with reprex v2.0.2</sup>
That one point crosses 180th meridian and apparently it silently transforms into EMPTY POINT
, which happens to be a valid geometry:
polygon_bbox_moll |> st_cast("POINT") |> st_transform(4326)
#> Simple feature collection with 5 features and 0 fields (with 1 geometry empty)
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 98.48587 ymin: -47.19364 xmax: 156.2191 ymax: -8.883331
#> Geodetic CRS: WGS 84
#> x
#> 1 POINT (124.3714 -47.19364)
#> 2 POINT EMPTY
#> 3 POINT (156.2191 -8.883331)
#> 4 POINT (98.48587 -8.883331)
#> 5 POINT (124.3714 -47.19364)
We can also plot Mollweide bbox with WGS84 graticule and transform 10x10 grid instead of 5-point polygon for better understanding. Transformed grid includes invalid polygons which plot()
can't handle, so we pass it through st_make_valid()
first:
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
polygon_bbox_moll <- st_bbox(c(xmin =9797297, xmax = 15540552, ymin = -5576498, ymax = -1096691), crs = st_crs("ESRI:54009")) |>
st_as_sfc() |>
st_as_sf()
plot(polygon_bbox_moll, graticule = TRUE, key.pos = NULL, axes = TRUE, las = 1)
polygon_bbox_moll |>
st_make_grid(n = 10) |>
st_transform(4326) |>
st_make_valid() |>
plot(graticule = TRUE, key.pos = NULL, axes = TRUE, las = 1)
Created on 2023-08-24 with reprex v2.0.2