Maps doesn't register weird shapes

Question

I'm working with one of my professors on some research aimed toward bettering the current methods of carbon accounting. We noticed that many of the locations for point sources were defaulted to the centroid of the county it was in (this is specific to the US at the moment, though it will be applied globally) if there was no data on the location.

So I'm using R to to address the uncertainty associated with these locations. My code takes the range of longitude and latitude for a county and plots 10,000 points. It then weeds out the points that are not in the county and take the average of the leftover points to locate the centroid. My goal is to ultimately take the difference between these points and the centroid to find the spacial uncertainty for point sources that were placed in the centroid.

However, I'm running into problems with coastal regions. My first problem is that the maps package ignores islands (the barrier islands for example) as well as other disjointed county shapes, so the centroid is not accurately reproduced when the points are averaged. My second problem is found specifically with Currituck county (North Carolina). Maps seems to recognize parts of the barrier islands contained in this county, but since it is not continuous, the entire function goes all wonky and produces a bunch of "NAs" and "Falses" that don't correspond with the actual border of the county at all.

(The data for the centroid is going to be used in other areas of the research which is why it's important we can accurately access all counties.)

Is there any way around the errors I'm running into? A different data set that could be read in, or anything of the sort? Your help would be greatly appreciated. Let me know if there are any questions about what I'm asking, and I'll be happy to clarify.

CODE:

ggplot2 helps

SOME TROUBLE COUNTIES: north carolina, currituck & massachusetts,dukes

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

First, given your description of the problem, I would probably invest a lot of effort to avoid this business of locations defaulting to the county centroids - that's the right way to solve this problem.

Second, if this is a research project, I would not use the built in maps in R. The USGS National Atlas website has excellent county maps of the US. Below is an example using Currituck County in NC.

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()

Finally, another way to do this (which I'd recommend, actually), is to just calculate the area weighted centroid. This is equivalent to what you are approximating, is more accurate, and much faster.

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

You'll find that as you increase the number of points in the points-weighted centroid the result gets closer to the area-weighted centroid.