rgisr-sf

Why does re-projecting a rectangle polygon result in a triangle using sf in R?


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>


Solution

  • 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