为什么在 R 中使用 sf 重新投影矩形多边形会产生三角形？

Question

我正在使用 R 包 sf。我在 Mollweide 投影中有一个澳大利亚专属经济区的边界框，我将其制作成 sf 多边形，然后重新投影到长纬度 (EPSG 4326) 坐标系。这会产生一个三角形多边形，如以下表示所示。我知道重新投影会扭曲多边形的矩形形状，但我不明白为什么它会产生三角形？

使用

st_coordinates

查看2个多边形的坐标，我可以看到经纬度中的多边形只有4个坐标，而Mollweide投影中的多边形有5个；但我不明白为什么删除了一个坐标。

任何帮助表示赞赏。

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

^{创建于 2023-08-24，使用 reprex v2.0.2}

Answer 1

那个点穿过 180 经度，显然它默默地转变为

EMPTY POINT

，这恰好是一个有效的几何图形：

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)

我们还可以使用 WGS84 标线绘制 Mollweide bbox，并变换 10x10 网格而不是 5 点多边形，以便更好地理解。变换后的网格包含

plot()

无法处理的无效多边形，因此我们首先将其传递给

st_make_valid()

：

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)

^{创建于 2023-08-24，使用 reprex v2.0.2}