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

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也省略NAs)来解决的问题。

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

...Instead提纲还是搞砸了，只是以一种不同的方式！

​

怎么了？

为了完整起见，下面是我用来掩盖国家边界以外高度数据的代码的其余部分：

# 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对其重新排序，则获得：

​

如果在使用函数SAml重新排序后使用基于sort_points的多边形，则可以获得：

​

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)