21

Rで緯度と経度の座標を州コードに変換する高速な方法はありますか? ルックアップ テーブルとして zipcode パッケージを使用してきましたが、緯度/経度の値を大量にクエリすると遅すぎます。

Rにない場合、Googleジオコーダーまたは他のタイプの高速クエリサービスを使用してこれを行う方法はありますか?

ありがとう！

4

5 に答える 5

45

ここには 2 つのオプションがあり、1 つはsfを使用し、もう 1 つはspパッケージ関数を使用します。sfは、空間データを分析するためのより近代的な (そして、ここでは 2020 年に推奨される) パッケージですが、それでもなお有用である場合に備えて、sp関連の関数でこれを行う方法を示す元の 2012 年の回答を残しています。

方法 1 (sf を使用):

library(sf)
library(spData)

## pointsDF: A data.frame whose first column contains longitudes and
##           whose second column contains latitudes.
##
## states:   An sf MULTIPOLYGON object with 50 states plus DC.
##
## name_col: Name of a column in `states` that supplies the states'
##           names.
lonlat_to_state <- function(pointsDF,
                            states = spData::us_states,
                            name_col = "NAME") {
    ## Convert points data.frame to an sf POINTS object
    pts <- st_as_sf(pointsDF, coords = 1:2, crs = 4326)

    ## Transform spatial data to some planar coordinate system
    ## (e.g. Web Mercator) as required for geometric operations
    states <- st_transform(states, crs = 3857)
    pts <- st_transform(pts, crs = 3857)

    ## Find names of state (if any) intersected by each point
    state_names <- states[[name_col]]
    ii <- as.integer(st_intersects(pts, states))
    state_names[ii]
}

## Test the function with points in Wisconsin, Oregon, and France
testPoints <- data.frame(x = c(-90, -120, 0), y = c(44, 44, 44))
lonlat_to_state(testPoints)
## [1] "Wisconsin" "Oregon"    NA

より高い解像度の州境界が必要な場合は、または他の手段sfを使用して、独自のベクター データをオブジェクトとして読み込みます。sf::st_read()便利なオプションの 1 つは、rnaturalearthパッケージをインストールし、それを使用してrnaturalearthhiresから状態ベクトル レイヤーを読み込むことです。次に、lonlat_to_state()ここに示すように定義したばかりの関数を使用します。

library(rnaturalearth)
us_states_ne <- ne_states(country = "United States of America",
                          returnclass = "sf")
lonlat_to_state(testPoints, states = us_states_ne, name_col = "name")
## [1] "Wisconsin" "Oregon"    NA

非常に正確な結果を得るには、GADM が管理する米国の行政境界を含むジオパッケージをこのページからダウンロードできます。次に、州の境界データを読み込み、次のように使用します。

USA_gadm <- st_read(dsn = "gadm36_USA.gpkg", layer = "gadm36_USA_1")
lonlat_to_state(testPoints, states = USA_gadm, name_col = "NAME_1")
## [1] "Wisconsin" "Oregon"    NA

方法 2 (sp を使用):

以下は、米国本土 48 州内の緯度と経度の data.frame を取得し、各ポイントについて、それが位置する州を返す関数です。

関数のほとんどは、パッケージ内の関数に必要な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

lonlat_to_state_sp <- 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(states_sp@polygons, function(x) x@ID)
    stateNames[indices]
}

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

lonlat_to_state_sp(testPoints)
[1] "wisconsin" "oregon" # IT WORKS
于 2012-01-06T00:39:24.740 に答える
3

?oversp パッケージを参照してください。

州の境界をSpatialPolygonsDataFrameとして持つ必要があります。

于 2012-01-06T00:02:28.933 に答える