r - map 不会记录奇怪的形状

Question

我正在与我的一位教授合作开展一些旨在改进当前碳核算方法的研究。我们注意到，如果没有关于位置的数据，许多点源的位置默认为它所在县的质心(目前这是美国特有的，但将在全局范围内应用)。

所以我使用 R 来解决与这些位置相关的不确定性。我的代码采用一个县的经度和纬度范围并绘制 10,000 个点。然后剔除不在县内的点，取剩余点的平均值来定位质心。我的目标是最终利用这些点与质心之间的差异来找到放置在质心中的点源的空间不确定性。

但是，我在沿海地区遇到了问题。我的第一个问题是 map 包忽略了岛屿(例如障碍岛)以及其他不连贯的县形状，因此在对点进行平均时无法准确再现质心。我的第二个问题是专门针对 Currituck 县(北卡罗来纳州)发现的。 map 似乎识别出该县的部分屏障岛，但由于它不是连续的，整个功能变得不稳定并产生一堆与实际边界不符的“NAs”和“Falses”一个县。

(质心的数据将用于研究的其他领域，这就是我们能够准确访问所有县的重要性。)

有什么方法可以解决我遇到的错误吗？可以读入的不同数据集，或任何类似的数据集？非常感谢您的帮助。如果对我的问题有任何疑问，请告诉我，我很乐意澄清。

代码:

ggplot2 帮助

一些麻烦县:北卡罗来纳州、柯里塔克和马萨诸塞州、杜克斯

library(ggplot2)
library(maps)     # package has maps
library(mapproj)  # projections
library(sp)
    
WC <- map_data('county','north carolina,currituck') #calling on county
   
p <- ggplot(data = WC, aes(x = long, y = lat)) #calling on latitude and longitude
   
p1 <- p + geom_polygon(fill = "lightgreen") + theme_bw() + 
  coord_map("polyconic") + coord_fixed() #+ labs(title = "Watauga County")
 
p1
   
### range for the long and lat
    
RLong <- range(WC$long)
RLong
   
RLat <- range(WC$lat)
RLat
    
### Add some random points
   
n <- 10000  
RpointsLong <- sample(seq(RLong[1], RLong[2], length = 100), n, replace = TRUE)
RpointsLat <- sample(seq(RLat[1], RLat[2], length = 100), n, replace = TRUE)
DF <- data.frame(RpointsLong, RpointsLat)
head(DF)
   
p2<-p1 + geom_point(data = DF, aes(x = RpointsLong, y = RpointsLat))
p2
  
#      Source:
# http://www.nceas.ucsb.edu/scicomp/usecases/GenerateConvexHullAndROIForPoints
    
inside <- map.where('county',RpointsLong,RpointsLat)=="north carolina,currituck"
inside[which(nchar(inside)==2)] <- FALSE
inside
   
g<-inside*DF
g1<-subset(g,g$RpointsLong!=0)
g1
   
CentLong<-mean(g1$RpointsLong)
CentLat<-mean(g1$RpointsLat)
Centroid<-data.frame(CentLong,CentLat)
Centroid

p1+geom_point(data=g1, aes(x=RpointsLong,y=RpointsLat)) #this maps all the points inside county
p1+geom_point(data=Centroid, aes(x=CentLong,y=CentLat))

Answer 1

首先，根据您对问题的描述，我可能会投入大量精力来避免这种位置业务默认为县质心 - 这是解决此问题的正确方法。

其次，如果这是一个研究项目，我不会使用 R 中的内置 map 。USGS National Atlas网站有优county maps of the US .以下是使用北卡罗来纳州柯里塔克县的示例。

library(ggplot2)
library(rgdal)    # for readOGR(...)
library(rgeos)    # for gIntersection(...)

setwd("< directory contining shapefiles >")
map   <- readOGR(dsn=".",layer="countyp010")
NC    <- map[map$COUNTY=="Currituck County" & !is.na(map$COUNTY),]
NC.df <- fortify(NC)

bbox  <- bbox(NC)
x   <- seq(bbox[1,1],bbox[1,2],length=50)   # longitude
y   <- seq(bbox[2,1],bbox[2,2],length=50)   # latitude
all <- SpatialPoints(expand.grid(x,y),proj4string=CRS(proj4string(NC)))
pts <- gIntersection(NC,all)                # points inside the polygons
pts <- data.frame(pts@coords)               # ggplot wants a data.frame
centroid <- data.frame(x=mean(pts$x),y=mean(pts$y))
ggplot(NC.df)+
  geom_path(aes(x=long,y=lat, group=group),    colour="grey50")+      
  geom_polygon(aes(x=long,y=lat, group=group), fill="lightgreen")+
  geom_point(data=pts, aes(x,y), colour="blue")+
  geom_point(data=centroid, aes(x,y), colour="red", size=5)+
  coord_fixed()

最后，另一种方法(实际上我会推荐)是只计算面积加权质心。这相当于您正在近似的值，更准确，也更快。

polys <- do.call(rbind,lapply(NC@polygons[[1]]@Polygons,
                              function(x)c(x@labpt,x@area)))
polys <- data.frame(polys)
colnames(polys) <- c("long","lat","area")
polys$area      <- with(polys,area/sum(area))
centr           <- with(polys,c(x=sum(long*area),y=sum(lat*area)))
centr       # area weighted centroid
#           x         y 
#   -76.01378  36.40105 
centroid    # point weighted centroid (start= 50 X 50 points)
#           x        y
# 1 -76.01056 36.39671

您会发现，随着点加权质心中点数的增加，结果会更接近面积加权质心。