我如何获得根据给定地点创建边界的纬度和经度？

Question

假设我的数据是这样的：

 daerah longitude   latitude    wilayah
    a   103.0059509 1.736281037   z
    a   103.0055008 1.736822963   z
    a   103.0049973 1.737220049   z
    a   103.0044479 1.737781048   z
    a   103.0041733 1.737781048   z
    a   103.003891  1.738060951   z
    a   103.0022202 1.738055944   z
    a   103.0019455 1.738332033   z
    a   103.0013885 1.738332033   z
    a   103.0011139 1.738610029   z
    a   103.0008316 1.738610029   z
    a   103.0005569 1.738891006   z
    a   103.000267  1.738891006   z
    a   103         1.738610029   z
    b   102.9966965 2.316540003   z
    b   102.9990997 2.315969944   z
    b   103.0125961 2.307929039   z
    b   103.0151978 2.306900978   z
    b   103.0171967 2.305169106   z
    b   103.0181961 2.30298996    z
    b   103.0189972 2.300110102   z
    b   103.0190964 2.29734993    z
    b   103.0169983 2.290781021   z
    b   102.9596024 2.197421074   z
    c   102.82444   2.365111113   z
    c   102.8239212 2.359646082   z
    c   102.8092346 2.338672638   z
    c   102.7966537 2.315601826   z
    c   102.7987518 2.290433407   z
    c   102.7987518 2.252681017   z
    c   102.7777786 2.225415468   z
    c   102.7421188 2.189760447   z
    c   102.7064667 2.183468342   z
    c   102.6708145 2.160397291   z
    c   102.6204758 2.137326479   z

daerah a,b,c 在wilayah z 中。
那么如何获得 wilayah 的经纬度呢？
经纬度应按照daerah点形成1个完整的形状。

任何帮助将不胜感激。提前谢谢你。

Answer 1

以下脚本中的data 是您的示例数据，其结构如下：

str(data)
# 'data.frame': 35 obs. of  4 variables:
#  $ daerah   : Factor w/ 3 levels "a","b","c": 1 1 1 1 1 1 1 1 1 1 ...
#  $ longitude: num  103 103 103 103 103 ...
#  $ latitude : num  1.74 1.74 1.74 1.74 1.74 ...
#  $ wilayah  : Factor w/ 1 level "z": 1 1 1 1 1 1 1 1 1 1 ...

使用包spatstat，您可以从所有数据创建一个多边形。

# install.packages("spatstat", dependencies = TRUE)
library(spatstat)

poly <- owin(xrange=range(data$longitude), 
     yrange=range(data$latitude), 
     poly=list(x=data$longitude, y=data$latitude))
plot(poly, axes=TRUE, las=1, col="lightblue")
points(data$longitude, data$latitude)

或者你可以更喜欢包sp，特别是如果点没有排序，以获得相同的结果（编辑：基于cmets添加）。

create.sp <- function(data, col.name, map.name){
    temp <- data[which(data[,col.name]==map.name),]
    poly <- Polygons(list(Polygon(temp[,c("longitude", "latitude")])), map.name)
    return(list(poly))
}

polys <- lapply(levels(data$wilayah), 
    FUN=function(x) SpatialPolygons(create.sp(data=data, col.name="wilayah", map.name=x)))

plot(polys[[1]], axes=TRUE, las=1, col="lightblue")

为每个 daerah 分别创建多边形并使用它们来表示一个 wilayah。 （编辑：基于 cmets 进行概括）。

library(sp)

# creates a single polygon from a character string of one daerah name
create.poly <- function(daerah){
    temp <- data[which(data$daerah==daerah),]
    return(Polygon(temp[,c("longitude", "latitude")]))
 }

# creates named polygons for all daerahs in the data 
# based on character string of their names
all.polys <- function(daerahs){
    polys <- list()
    for(i in 1:length(daerahs)){
       daerah <- levels(data$daerah)[i]
       polys[[i]] <- Polygons(list(create.poly(daerah)), daerah)
    }
    return(polys)
}

polys <- SpatialPolygons(all.polys(levels(data$daerah)), 1:length(levels(data$daerah)))

plot(polys, axes=TRUE, las=1, col="lightblue")

或者您可以使用包alphahull 从您的数据点绘制一个凹壳。

library(alphahull)

phull <- ahull(data[,c("longitude","latitude")], alpha = 1)
pshap <- ashape(data[,c("longitude","latitude")], alpha = .3)

par(mfrow=c(1,2))
plot(phull, main="alpha hull")
plot(pshap, main="alpha shape")

凹壳不代表任何类型的spatial 对象，但这些函数对于在一组点周围绘制区域非常方便。

您的问题的理想解决方案可能是ashape 边界，因为它不会像早期的选项那样与自身相交。不过不知道怎么解决this problem。