我正在与我的一位教授一起进行一些旨在改善当前碳会计方法的研究。我们注意到，如果没有点数据源，则许多点源位置默认为所在县的质心(目前特定于美国，尽管将在全球范围内应用)。
因此，我使用R来解决与这些位置相关的不确定性。我的代码接受一个县的经度和纬度范围，并绘制10,000点。然后，它清除掉不在县中的点，并取剩余点的平均值来定位质心。我的目标是最终求出这些点和质心之间的差，以找到放置在质心中的点源的空间不确定性。
但是，我遇到了沿海地区的问题。我的第一个问题是，地图包忽略了岛(例如，障碍岛)以及其他不连贯的县形状，因此，在对点进行平均时，质心无法得到精确再现。我的第二个问题是在Currituck县(北卡罗来纳州)发现的。地图似乎可以识别出该县包含的部分障碍岛，但由于它不是连续的，因此整个功能全都变得很奇怪，并产生了与实际边界不符的“NA”和“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))
最佳答案
首先，考虑到您对问题的描述，我可能会花费很多精力来避免这种位置默认为县质心的问题-这是解决此问题的正确方法。
其次，如果这是一个研究项目，我不会使用R中的内置地图。USGSNational Atlas网站的county maps of the US非常出色。以下是在北卡罗来纳州使用Currituck County的示例。
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
您会发现，随着增加点加权质心中的点数，结果将更接近面积加权质心。