はじめに
はじめまして！データサイエンティストの渡邊です。株式会社GEOTRA(@GEOTRA )にてGIS（Geographic Information System：地理情報システム）データを用いた分析を行っています。
本記事ではgeopandasで対象のPOLYGONに含まれるLINESTRINGを抽出する方法を紹介します。
背景
- 道路に関するGISデータを扱う際、与えられるデータはしばしば膨大になり、分析対象の絞り込みが必要になることがあります
- 例えば、分析対象のエリアが決まっている場合はそのエリアを通る道路のみを抽出することで計算量の削減に繋げることができます
やること
- 今回は道路データのサンプルデータとして国土数値情報ダウンロードサイトの道路データを用いたPythonの分析例を紹介します
- 用いるのは、道路法に基づく高速自動車国道、一般国道、都道府県道、市町村道及び特例都道等、全国の道路について、位置（線）、路線名等を整備したデータとなります(以下サンプル)
| N01_001 | N01_002 | N01_003 | N01_004 | geometry |
|---|---|---|---|---|
| 3 | 東京所沢線 | None | None | LINESTRING (139.4967777800000022 35.7993055600... |
| 3 | 東京所沢線 | None | None | LINESTRING (139.4920555600000114 35.7984722200... |
| 3 | 東京所沢線 | None | None | LINESTRING (139.4901111099999866 35.7984722200... |
| 3 | 東京所沢線 | None | None | LINESTRING (139.4876111099999889 35.7979166699... |
| 3 | 東京所沢線 | None | None | LINESTRING (139.4695555600000034 35.7943055599... |
- 道路を表す
LINESTRINGと対象エリアを表すPOLYGONについて以下の場合を考えます- 対象エリアと少しでも被る道路を取得する場合
- 対象エリアに完全に含まれる道路を取得する場合
- 今回は東京都千代田区を対象エリアとします
- ポリゴンの作り方は以下参照
対象エリアと少しでも被る道路を取得する場合
-
geopandas.geometry.intersectsを用います- 交差する場合は
True、そうでなければFalseが返されます
- 交差する場合は
サンプルコード
from keplergl import KeplerGl # 0.3.2
import geopandas as gpd # 0.10.2
from shapely.geometry import Polygon # 1.8.0
# 東京都の道路データ(世界測地系) https://nlftp.mlit.go.jp/ksj/gmlold/datalist/gmlold_KsjTmplt-N01.html
gdf_road = gpd.read_file('N01-07L-13-01.0a_GML/N01-07L-2K-13_Road.shp')
# 千代田区周辺のPOLYGON
polygon = Polygon([[139.74598303331285,35.70292105511908],[139.7312914052023,35.689224275705776],[139.73001713133573,35.685266990373435],[139.7335401237903,35.67875226236136],[139.73691320167327,35.67856959876946],[139.73931183483475,35.67357663202999],[139.74793192275757,35.66955767581042],[139.749655940341,35.67077556270092],[139.75497790766673,35.66967946533563],[139.76427261116578,35.67241968052733],[139.770494065927,35.67808249381528],[139.77109372421742,35.685145224635534],[139.76989440763668,35.6870325726985],[139.77064398050055,35.689346035218314],[139.78278706087795,35.694825020913306],[139.78241227444698,35.69786873921977],[139.7732674855196,35.70486885044058],[139.77004432221017,35.70529492432322],[139.7663714151821,35.70182539933901],[139.76824534733908,35.69975578622787],[139.7620988498636,35.70109495378689],[139.74598303331285,35.70292105511908]])
gdf_polygon = gpd.GeoDataFrame(geometry=[polygon])
# 対象エリアと少しでも被る道路を取得
# 各LINESTRINGが該当POLYGONと交差するか否かをbool値で返す
gdf_road['intersection'] = gdf_road['geometry'].intersects(polygon)
# keplerglで可視化
kepler_map = KeplerGl(height=800)
kepler_map.add_data(data=gdf_road, name='Road')
kepler_map.add_data(data=gdf_polygon, name='千代田区')
kepler_map
intersectionの値によって色分けするとこのようになります
対象エリアに完全に含まれる道路を取得する場合
-
geopandas.geometry.withinまたはshapely.Polygon.containsを用います-
POLYGON内に含まれる場合はTrueそうでなければFalseが返されます
-
サンプルコード
from keplergl import KeplerGl # 0.3.2
import geopandas as gpd # 0.10.2
from shapely.geometry import Polygon # 1.8.0
# 東京都の道路データ(世界測地系) https://nlftp.mlit.go.jp/ksj/gmlold/datalist/gmlold_KsjTmplt-N01.html
gdf_road = gpd.read_file('N01-07L-13-01.0a_GML/N01-07L-2K-13_Road.shp')
# 千代田区周辺のPOLYGON
polygon = Polygon([[139.74598303331285,35.70292105511908],[139.7312914052023,35.689224275705776],[139.73001713133573,35.685266990373435],[139.7335401237903,35.67875226236136],[139.73691320167327,35.67856959876946],[139.73931183483475,35.67357663202999],[139.74793192275757,35.66955767581042],[139.749655940341,35.67077556270092],[139.75497790766673,35.66967946533563],[139.76427261116578,35.67241968052733],[139.770494065927,35.67808249381528],[139.77109372421742,35.685145224635534],[139.76989440763668,35.6870325726985],[139.77064398050055,35.689346035218314],[139.78278706087795,35.694825020913306],[139.78241227444698,35.69786873921977],[139.7732674855196,35.70486885044058],[139.77004432221017,35.70529492432322],[139.7663714151821,35.70182539933901],[139.76824534733908,35.69975578622787],[139.7620988498636,35.70109495378689],[139.74598303331285,35.70292105511908]])
gdf_polygon = gpd.GeoDataFrame(geometry=[polygon])
# ========== 書き換え箇所 ==========
# 対象エリアに完全に含まれる道路を取得
# 各LINESTRINGが該当POLYGON内にあるか否かをbool値で返す
gdf_road['within'] = gdf_road['geometry'].within(polygon)
gdf_road['contains'] = [polygon.contains(g) for g in gdf_road['geometry']]
# keplerglで可視化
kepler_map = KeplerGl(height=800)
kepler_map.add_data(data=gdf_road, name='Road')
kepler_map.add_data(data=gdf_polygon, name='千代田区')
kepler_map
within/contains の値によって色分けするとこのようになります
おわりに
最後までお読みくださりありがとうございました。
GEOTRAのQiitaでは、引き続き地理空間情報に関する情報を中心に技術系のテクニックを発信していきます。今後も皆さんのお役に立てるコンテンツを配信できればと思っておりますので、皆様のいいね・フォローをお待ちしています！
また、こちらの記事を読んでGEOTRAに興味を持ってくださった方は、是非Wantedlyから会社紹介記事もご覧ください！GEOTRAは、学生インターンを含め、現在一緒に働く仲間を募集しております。興味をお持ちの方は、Wantedly上の募集ページやストーリー、各種サイトをご確認ください！
企業HP： https://www.geotra.jp/
note： https://note.com/2022geotra/
Wantedly： https://www.wantedly.com/stories/s/geotra_introduction
LinkedIn： https://www.linkedin.com/company/%E6%A0%AA%E5%BC%8F%E4%BC%9A%E7%A4%BEgeotra/