对坐标进行排序以创建多边形会产生混乱的结果

Question

如何从一组坐标开始获得具有南非形状的多边形，这些坐标表示在 data.frame 中组织的国家边界？我已经很接近了，但 我无法创建具有正确国家形状的多边形。

我正在尝试绘制类似于以下内容的南非海拔地图：

但是在国家边界之外的高度信息被掩盖了。这是我一直在使用的代码：

# Load relevant packages:
library(elevatr)
library(raster)
library(sf)
library(sp)
library(fasterize)
library(maps)
library(mapdata)

SA_sf <- map("worldHires", "South Africa")
str(SA_sf); head(SA_sf)
map(SA_sf) # OK
# Extract coordinates
SA_coords <- data.frame(x=SA_sf$x, y=SA_sf$y)
plot(SA_coords, type="l") # just as it should be

在这里，我只是简单地进行子集化以仅保留连续的土地和莱索托：

SAml <- subset(SA_coords, x>16 & x<35 & y > -40 & y < -20)
plot(SAml, type="l")

大纲搞砸了：

以下应该通过按顺时针顺序对坐标进行排序来解决问题（sort_points 也省略了 NA）...

# library(devtools)
# install_github("skgrange/gissr")
library(gissr)
SAml <- sort_points(SAml, y = "y", x = "x", clockwise = T)
plot(SAml, type="l")

...相反，轮廓仍然混乱，只是方式不同！

怎么了？

为了完整起见，下面是我用来屏蔽国家边界以外的海拔数据的其余代码：

# Create a SpatialPoints object
prj_dd <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
examp_sp <- SpatialPoints(SAml, proj4string = CRS(prj_dd))

# get raster file of altitude:
sp_elev_epqs <- get_elev_raster(examp_sp, z=9)
# this takes a while to download, so once it's done it makes sense to save it on your computer:
# writeRaster(sp_elev_epqs, filename="~/Desktop/SA_mainland.grd")

# first, close your polygon 
# (first and last points must be identical)
SAml <- rbind(SAml, SAml[1,])
# Use various functions from the sf package to do the magic:
poly <- st_sf(st_sfc(st_polygon(list(as.matrix(SAml)))), crs = 4326)
str(poly)
# library(fasterize)
poly$POLYID <- 1:nrow(poly)
elevation_data <- sp_elev_epqs
polymap <- fasterize(poly, elevation_data, field = "POLYID")
## mask out elevation where there's no polygon
elevation_data[is.na(values(polymap))] <- NA
# or use:
# elevation_data <- mask(x=sp_elev_epqs, mask= poly)
plot(elevation_data, col = gray.colors(80, start = 0.9, end = 0.1, gamma = 2.2, alpha = NULL))
lines(SA_coords)

如果我使用基于SAml 的多边形而不使用函数sort_points 对其重新排序，我会得到：

如果我在使用函数sort_points 重新排序后使用基于SAml 的多边形，我会得到：

Answer 1

我会先创建一个 SpatialPolygons 对象

library(raster)
library(mapdata)
require(maptools)

x <- map("worldHires", "South Africa", fill=TRUE)
y <- map2SpatialPolygons(x, IDs=1:length(x$names), proj4string=CRS("+proj=longlat +datum=WGS84"))

在这种情况下，您可以简单地保留最大的多边形

a <- area(y)
b <- y[which.max(a)]
plot(b)

如果你想要一个 data.frame

g <- geom(b)

一个更简单的替代方法是

library(raster)
sa <- getData("GADM", country="South Africa", level=0)