如何在R中无问题地将经纬度转换为FIPS县代码？

Question

我有成千上万的经纬度点，我用这段代码把它们转换成FIPS县代码。

latlong2county <- function(pointsDF) {
# Prepare SpatialPolygons object with one SpatialPolygon
# per state (plus DC, minus HI & AK)
states <- map('county', fill=TRUE, col="transparent", plot=FALSE)
IDs <- sapply(strsplit(states$names, ":"), function(x) x[1])
states_sp <- map2SpatialPolygons(states, IDs=IDs,
                               proj4string=CRS("+proj=longlat +datum=wgs84"))

# Convert pointsDF to a SpatialPoints object 
pointsSP <- SpatialPoints(pointsDF, 
                        proj4string=CRS("+proj=longlat +datum=wgs84"))

# Use 'over' to get _indices_ of the Polygons object containing each point 
indices <- over(pointsSP, states_sp)

# Return the state names of the Polygons object containing each point
stateNames <- sapply(states_sp@polygons, function(x) x@ID)
stateNames[indices]
}

data(county.fips)
latlon <- data.frame(all$lon,all$lat)
county<-latlong2county(latlon)
fips<-with(county.fips, fips[match(county, polyname)])

但是，这个方法的问题是并不是所有的经纬度都能转换成FIPS代码。有些合法的县居然得到了NA值。有没有人能帮帮我？

或者简单地推荐一个不同的解决方案？

Answer 1

虽然这个解决方案不是基于Python的，但你可以使用Ruby和geokit这个库来提前处理你的数据，然后通过csv文件或其他方式把它传给R。

想要了解更多的步骤，可以查看这个文章。