2012-01-05 2 views
15

Есть ли быстрый способ преобразования координат широты и долготы в коды состояний в R? Я использую пакет zipcode в качестве таблицы поиска, но он слишком медленный, когда я запрашиваю множество значений lat/long.Широта долготы Координаты с кодом состояния в R

Если нет в R, то есть какой-либо способ сделать это, используя геокодер google или любой другой тип службы быстрого запроса?

Спасибо!

+0

см также мой ответ здесь, используя 'ggmap :: revgeocode': https://stackoverflow.com/questions/46150851/how-to-get-california- уездное местоположение от однолатового itude-and-longitude-information/46151310 # 46151310 –

ответ

32

Вот функция, которая берет data.frame lat-longs в нижних 48 состояниях и для каждой точки возвращает состояние, в котором оно расположено.

Большая часть функции просто готовит SpatialPoints и SpatialPolygons объектов, необходимых over() функции в sp пакете, который делает реальные тяжелый расчет «пересечения» точек и полигонов:

library(sp) 
library(maps) 
library(maptools) 

# The single argument to this function, pointsDF, is a data.frame in which: 
# - column 1 contains the longitude in degrees (negative in the US) 
# - column 2 contains the latitude in degrees 

latlong2state <- function(pointsDF) { 
    # Prepare SpatialPolygons object with one SpatialPolygon 
    # per state (plus DC, minus HI & AK) 
    states <- map('state', 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([email protected], function(x) [email protected]) 
    stateNames[indices] 
} 

# Test the function using points in Wisconsin and Oregon. 
testPoints <- data.frame(x = c(-90, -120), y = c(44, 44)) 

latlong2state(testPoints) 
[1] "wisconsin" "oregon" # IT WORKS 
+2

Мне пришлось изменить wgs84 на WGS84, чтобы этот пример работал. – lever

+0

@lever Спасибо, что указали это. Интересно, когда (и где) было изменено. В любом случае, я теперь отредактировал, чтобы исправить это. –

2

Подробнее см. В пакете sp. Вы должны иметь границы состояний как SpatialPolygonDataFrame.

4

You может сделать это в несколько строк R.

library(sp) 
library(rgdal) 
#lat and long 
Lat <- 57.25 
Lon <- -9.41 
#make a data frame 
coords <- as.data.frame(cbind(Lon,Lat)) 
#and into Spatial 
points <- SpatialPoints(coords) 
#SpatialPolygonDataFrame - I'm using a shapefile of UK counties 
counties <- readOGR(".", "uk_counties") 
#assume same proj as shapefile! 
proj4string(points) <- proj4string(counties) 
#get county polygon point is in 
result <- as.character(over(points, counties)$County_Name) 
+0

Спасибо! Это было еще проще :) –