データフレーム上の緯度と経度を空間オブジェクトに変換する-🦉sfと🐍geopandasの例-
地理空間データを取り扱う際は、はじめにgeojsonやshapeファイルで受け取ることが一般的かと思います。 あるいはポイントデータの場合には、緯度と経度の値が各列に記録されるcsvなどの表形式のテキストファイルを起点とすることもあります。
前者のような地理空間データであれば、 R等のアプリケーションで読み込めば自動的に地理空間データとみなしてくれますが、 後者の場合はそうではありません。 座標の値はあくまでも数値です。なのでこうしたデータを地理空間データとして扱えるようにするには変換作業が必要となります。
今回の記事では、RおよびPythonでのデータフレームに記録された緯度経度の列を変換し、空間オブジェクトとして扱えるようにする方法を紹介します。 空間オブジェクトの形式としてRではsf、Pythonではgeopandasを対象にします。地理空間データに変換しておくと、データの空間的な配置を可視化可能になったり空間演算が可能になり、分析の幅を広げることが期待できます。
対象データ
csvファイルに緯度と経度の値が格納されているデータの例として、国土交通省 位置参照情報ダウンロードサービス 提供の平成30年整備 茨城県つくば市の大字・町丁目レベルのデータを利用します。
ファイルをダウンロードしたら、まずはデータフレームとしてcsvファイルを読み込みましょう。
library(sf) library(readr)
df_isj08220 <- read_csv("08220_2018.csv", locale = locale(encoding = "cp932"), col_types = "___c_cdd__") %>% purrr::set_names(c("city", "street_lv1", "latitude", "longitude"))
データの一部を表示します。latitude, longitudeの列がそれぞれ独立しています。
city | street_lv1 | latitude | longitude |
---|---|---|---|
つくば市 | 赤塚 | 36.04268 | 140.1232 |
つくば市 | 明石 | 36.18207 | 140.0537 |
つくば市 | 旭 | 36.11289 | 140.0725 |
つくば市 | 安食 | 36.16775 | 140.0163 |
つくば市 | あしび野 | 35.96060 | 140.1122 |
つくば市 | 吾妻一丁目 | 36.08153 | 140.1129 |
sfパッケージ
sf
パッケージでデータフレームの座標情報を参照してポイントデータに変換するにはst_as_sf()
を使います。
第一引数に座標列を含む対象のデータフレーム、coords引数の値にX,Y座標となる列名を指定します。この時、座標参照系が判明している際はcrs引数によって明示しておきましょう。
sf_isj08220 <- df_isj08220 %>% sf::st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
sf_isj08220 ## Simple feature collection with 299 features and 2 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 140.0102 ymin: 35.95047 xmax: 140.1676 ymax: 36.22462 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## # A tibble: 299 x 3 ## city street_lv1 geometry ## <chr> <chr> <POINT [°]> ## 1 つくば市 赤塚 (140.1232 36.04268) ## 2 つくば市 明石 (140.0537 36.18207) ## 3 つくば市 旭 (140.0725 36.11289) ## 4 つくば市 安食 (140.0163 36.16775) ## 5 つくば市 あしび野 (140.1122 35.9606) ## 6 つくば市 吾妻一丁目 (140.1129 36.08153) ## 7 つくば市 吾妻二丁目 (140.1111 36.08443) ## 8 つくば市 吾妻三丁目 (140.1109 36.0889) ## 9 つくば市 吾妻四丁目 (140.1183 36.08722) ## 10 つくば市 天久保一丁目 (140.1087 36.0919) ## # … with 289 more rows
緯度経度の列を元にポイントデータが作成できました。地図上にマッピングしてみます。
mapview::mapview(sf_isj08220)
sfheadersパッケージ
続いてsfheaders
パッケージを使う方法です。sf::st_as_sf()
の実行結果と同じく、sfオブジェクトを作成します。
今回は対象がポイントデータなので、sf_point()
のx,y引数にそれぞれ経度と緯度の列名を与えます。
この時、keep引数の値をTRUEにしておくとデータフレームの他の列が残ります。FALSEではポイントのgeometry列のみしか残らないため注意です。
また、座標参照系もこのタイミングでは指定できないため、sf::st_set_crs()
で定義する必要が生じます。
df_isj08220 %>% sfheaders::sf_point(x = "longitude", y = "latitude", keep = TRUE) %>% st_set_crs(4326) %>% tibble::new_tibble(subclass = "sf", nrow = nrow(.))
sf_isj08220sfh ## Simple feature collection with 299 features and 2 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: 140.0102 ymin: 35.95047 xmax: 140.1676 ymax: 36.22462 ## z_range: zmin: NA zmax: NA ## m_range: mmin: NA mmax: NA ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## # A tibble: 299 x 3 ## city street_lv1 geometry ## <chr> <chr> <POINT [°]> ## 1 つくば市 赤塚 (140.1232 36.04268) ## 2 つくば市 明石 (140.0537 36.18207) ## 3 つくば市 旭 (140.0725 36.11289) ## 4 つくば市 安食 (140.0163 36.16775) ## 5 つくば市 あしび野 (140.1122 35.9606) ## 6 つくば市 吾妻一丁目 (140.1129 36.08153) ## 7 つくば市 吾妻二丁目 (140.1111 36.08443) ## 8 つくば市 吾妻三丁目 (140.1109 36.0889) ## 9 つくば市 吾妻四丁目 (140.1183 36.08722) ## 10 つくば市 天久保一丁目 (140.1087 36.0919) ## # … with 289 more rows
おまけ: sfからdfに戻す
データフレーム(df)の緯度経度情報をもとにsfオブジェクトを作成する例を述べましたが、反対にsfからデータフレームへ、地理空間情報を取り除く方法も紹介します。
作成したsfオブジェクトに対してsf::st_drop_geometry()
ないしsf::st_set_geometry(value = NULL)
でgeometry情報を除外することでdfに変換されます。
sf_isj08220 %>% st_drop_geometry() sf_isj08220 %>% st_set_geometry(value = NULL)
この時、ポイントの情報も失われるので、データフレームの列として残したい場合は予め列を追加する必要が生じます。
sf_isj08220 %>% dplyr::mutate(coords = purrr::map(geometry, ~ as.data.frame(sf::st_coordinates(.x)))) %>% tidyr::unnest(cols = coords) %>% st_drop_geometry() ## # A tibble: 299 x 4 ## city street_lv1 X Y ## * <chr> <chr> <dbl> <dbl> ## 1 つくば市 赤塚 140. 36.0 ## 2 つくば市 明石 140. 36.2 ## 3 つくば市 旭 140. 36.1 ## 4 つくば市 安食 140. 36.2 ## 5 つくば市 あしび野 140. 36.0 ## 6 つくば市 吾妻一丁目 140. 36.1 ## 7 つくば市 吾妻二丁目 140. 36.1 ## 8 つくば市 吾妻三丁目 140. 36.1 ## 9 つくば市 吾妻四丁目 140. 36.1 ## 10 つくば市 天久保一丁目 140. 36.1 ## # … with 289 more rows
Python
import pandas as pd import geopandas
df = pd.read_csv("08220_2018.csv", encoding = "cp932", usecols = [3, 5, 6, 7], dtype = {'市区町村名': str, '大字町丁目名': str, '緯度': float, '経度': float}) df = df.rename(columns = {'市区町村名': 'city', '大字町丁目名': 'street_lv1', '緯度': 'latitude', '経度': 'longitude'})
geopandas
gdf = geopandas.GeoDataFrame(df, geometry = geopandas.points_from_xy(df.longitude, df.latitude))
gdf.head(10)
city | street_lv1 | longitude | latitude | geometry |
---|---|---|---|---|
つくば市 | 赤塚 | 140.1 | 36.04 | POINT (140.12315 36.042677) |
つくば市 | 明石 | 140.1 | 36.18 | POINT (140.053737 36.182065) |
つくば市 | 旭 | 140.1 | 36.11 | POINT (140.072495 36.11289) |
つくば市 | 安食 | 140.0 | 36.17 | POINT (140.016345 36.167751) |
つくば市 | あしび野 | 140.1 | 35.96 | POINT (140.112185 35.960604) |
つくば市 | 吾妻一丁目 | 140.1 | 36.08 | POINT (140.112938 36.081534) |
つくば市 | 吾妻二丁目 | 140.1 | 36.08 | POINT (140.111144 36.084425) |
つくば市 | 吾妻三丁目 | 140.1 | 36.09 | POINT (140.11092 36.088898) |
つくば市 | 吾妻四丁目 | 140.1 | 36.09 | POINT (140.118261 36.087219) |
つくば市 | 天久保一丁目 | 140.1 | 36.09 | POINT (140.108708 36.091901) |
せっかくなのでこちらもマッピング。
import geopatra m = gdf.folium.plot(zoom = 10) m
いいですね!Enjoy!