国土数値情報道路データで道路の幅員をプロットしてみる

はじめに

なんとなく坂戸市の自宅の前の道路について知りたくなったので、道路の幅員情報をプロットするアプリをつくってみた。

幅員は色で識別される

• GitHubリポジトリ
  https://github.com/tatamyiwathy/SakadoRoads

道路データについて

道路データは国土数値情報ダウンロードサイトの「道路データ」からダウンロードした。

このページではダウンロードするファイルをメッシュで選択する。
しかし、私が欲しいのは市区町村ごとの道路のデータであるが、ダウンロードページでは市区町村の認識が困難である。
また、市区町村がメッシュで分割されてしまっていると、複数の道路データをダウンロードする必要がある。
そこで、市区町村を指定してそこにある道路だけを抽出してGeoJSON形式で出力するPythonスクリプトを作成した。

このスクリプトはPythonで書かれている。geopandasが必要なのであらかじめインストールしておく。

$ pip install geopandas

使い方

$ python create-road-dataset.py

引数	説明
boundary	市区町村名
output	出力ファイル名（GeoJSON形式）
-p	都道府県名（オプション）

第1引数で抽出したい市区町村、第2引数で出力ファイルを指定する。

第1引数で指定した市区町村に同一名称の別の市区町村が存在していると、それらすべてが抽出されてしまうので-pオプションで都道府県を指定する。

例：東京都府中市、広島県府中市

$ python create-road-dataset.py -p 広島県 府中市 road_dataset.geojson

広島県府中市だけが抽出される。
なお市区町村の指定が「府中」だけだと抽出に失敗する。「府中市」のように市区町村を明示する。

完了すると、実行した場所に第2引数で指定したファイルが作成される。
市区町村がメッシュで分割されていたとしても複数の道路データをマージしてgeojson形式のファイルを出力してくれる。

スクリプトの動作について

1. 市区町村の境界を得るため、国土数値情報ダウンロードサイト（行政区域データ）から「全国」のデータをダウンロードして、zip展開する
2. もし-pオプションがあるなら都道府県を指定して、行政区域を抽出する
3. boundaryで行政区域データを抽出する
4. total_boundsで市区町村のバウンダリーボックスを取得する
5. バウンダリーボックスで得られた経緯度を1次メッシュコードに変換する（市区町村がメッシュで分割されているとメッシュコードが複数できる）

   def get_mesh_code(latitude, longitude):
     """
     Calculate the mesh code for a given latitude and longitude.
   
     Args:
         latitude (float): Latitude in degrees.
         longitude (float): Longitude in degrees.
   
     Returns:
         str: Mesh code in the format "LLLLOO", where LLLL is the latitude code and OO is the longitude code.
     """
     lat_code = int(latitude * 1.5)
     lon_code = int(longitude) - 100
     return f"{lat_code:02d}{lon_code:02d}"
   

6. 5.で得られたメッシュコードを元に国土数値情報ダウンロードサイト（道路データ）からダウンロードしてzip展開する
7. 道路データが複数ならひとつに結合する
8. 4.で得られたバウンダリーボックスの外にある道路を除外する
9. 2.で得られた行政区域の境界線の外にある道路を除外する
   is_road_in_boundary()が道路を構成する経緯度毎に呼び出され判定を行う

   def is_point_in_boundary(point, boundary):
       """
       Check if a point is within a boundary.
   
       Args:
           point (shapely.geometry.Point): The point to check.
           boundary (geopandas.GeoDataFrame): The GeoDataFrame of the boundary.
   
       Returns:
           bool: True if the point is within the boundary, False otherwise.
       """
       # Check if the point is within the boundary's geometry
       return boundary.geometry.contains(point).any()
   
   
   def is_road_in_boundary(geometry, boundary):
       """
       Check if a road (LineString or MultiLineString) is within a boundary.
   
       Args:
           geometry (shapely.geometry.LineString or     shapely.geometry.MultiLineString): The geometry of the road.
   
       Returns:
           bool: True if any part of the road is within the boundary, False     otherwise.
       """
   
       if isinstance(geometry, LineString):
           lines = [geometry]
       elif isinstance(geometry, MultiLineString):
           lines = geometry.geoms
       else:
           return False
   
       return any(
           is_point_in_boundary(Point(longitude, latitude), boundary)
           for line in lines
           for longitude, latitude in line.coords
       )
   
   

10. outputで指示されたファイルにGeoJSON形式で出力する