Python Remind Memo - World Heritage Plotting in Python and ChatGPT

日付:

国土数値情報の地域資源・観光という項目に「世界文化遺産」、「世界自然遺産」というコンテンツがあったので、図化を試みてみました。(URL:https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A34-v1_2.html)

まずは、GIS情報（Shapeファイル）のダウンロードから始めますが、今回は手動でダウンロード→解凍→GoogleDriveにアップロードするのではなく、GoogleClabの中でこれらの手順を済ませてしまいたいと思いました。
調べてみると、GeoPandasにZIPファイルを直接読み込む機能があるようなので、以下のコードを実行しました。
!pip install geopandas
import geopandas
import pathlib
NOTEBOOK_PATH = pathlib.Path().resolve()
DATA_DIRECTORY = NOTEBOOK_PATH / "data"

print("####", DATA_DIRECTORY)

world_heritages = geopandas.read_file('https://nlftp.mlit.go.jp/ksj/gml/data/A34/A34-17/A34-180316_GML.zip')
world_heritages.head()


ちなみに、'https://nlftp.mlit.go.jp/ksj/gml/data/A34/A34-17/A34-180316_GML.zip'は国土数値情報の該当ページのHTMLコードを解読して、ZIPファイルの保存場所を調べました。
しかし、この結果は次のようなErrorになりました。
ERROR:fiona._env:`/vsizip/vsicurl/https://nlftp.mlit.go.jp/ksj/gml/data/A34/A34-17/A34-180316_GML.zip' not recognized as a supported file format.

ChatGPT先生にこのエラーの理由をお尋ねすると、以下のご教示をいただきました。

このエラーは、geopandas.read_file()メソッドが指定されたURLから直接ZIPファイルを読み込むことができないことを示しています。read_file()メソッドは、URLから直接ファイルを読み込む機能を持っていないため、このエラーが発生しました。

ZIPファイルをダウンロードして解凍するためには、先ほど紹介した方法を使用して手動で処理する必要があります。以下に、改めてZIPファイルをダウンロードして解凍する手順を示します。

python
import urllib.request
import zipfile

zip_url = "https://nlftp.mlit.go.jp/ksj/gml/data/A34/A34-17/A34-180316_GML.zip"

# ZIPファイルをダウンロード
urllib.request.urlretrieve(zip_url, "example.zip")

# ZIPファイルを解凍
with zipfile.ZipFile("example.zip", "r") as zip_ref:
    zip_ref.extractall("destination_folder")

上記のコードの最後の行の"destination_folder"の部分を、ZIPファイルを解凍したいディレクトリのパスに変更してください。これにより、ZIPファイルがダウンロードされ、指定したディレクトリに解凍されます。

もし、解凍後のデータを読み込むためにgeopandasを使用したい場合は、解凍されたファイルのパスを指定してgeopandas.read_file()を呼び出すことができます。

python
import geopandas

shapefile_path = "destination_folder/A34-180316_GML/A34-180316.shp"
world_heritages = geopandas.read_file(shapefile_path)
world_heritages.head()

上記のコードでは、解凍されたデータの中にある.shpファイルを指定してgeopandas.read_file()を使用しています。もちろん、解凍されたデータのファイル形式に応じて、適切なファイルを指定してください。


ここでありがちなエラー。
Error:unsupported operand type(s) for +: 'PosixPath' and 'str'
というのが発生しました。こちらもChatGPT先生にご質問したところ、


おっしゃるエラーは、PosixPath型（またはPath型）のオブジェクトと文字列型のオブジェクトを+演算子で結合しようとした際に発生するものです。Pathオブジェクトは、ファイルシステムのパスを操作するためのクラスであり、文字列とは異なる型です。

このエラーを解決するためには、Pathオブジェクトを文字列に変換してから結合する必要があります。以下のように、str()関数を使用してPathオブジェクトを文字列に変換してください。

python
import pathlib

# DATA_DIRECTORYがPosixPathオブジェクトとして定義されている場合
destination_folder = str(DATA_DIRECTORY) + "/A34-180316_GML"
shapefile_path = str(destination_folder) + "/A34-180316.shp"

または、Pathオブジェクトを/演算子で連結する方法もあります。

python
import pathlib

# DATA_DIRECTORYがPosixPathオブジェクトとして定義されている場合
destination_folder = DATA_DIRECTORY / "A34-180316_GML"
shapefile_path = destination_folder / "A34-180316.shp"

いずれの方法でも、Pathオブジェクトを文字列に変換するか、/演算子を使ってPathオブジェクト同士を結合することで、エラーを回避できます。


だそうです。Pathlibで作成したPathhはPosixPath型という特殊なTypeらしく、PosixPath型（またはPath型）のオブジェクトと文字列型のオブジェクトを+演算子で結合できないとのことで、一度PosixPathをstr()で文字列に変換する作業がいるようです。そんな面倒な手続きを踏まずにPathlibが直接文字列を吐き出せばいいような気がしますが、何かりゆうがあるのでしょうね？
ちょっとChatGPT先生に聞いてみました。


pathlibライブラリは、様々な操作を行うためにパスを表現するための高水準なAPIを提供します。その際に、PosixPath型やWindowsPath型など、環境に応じた具体的なパスオブジェクトの型を返すことがあります。

pathlibの設計哲学の一つは、環境に依存しないパス表現を提供することです。つまり、pathlibはWindowsやUnix系の異なる環境で使用されることを想定しており、それぞれの環境に適したパス表現を返すことができます。例えば、PosixPathはUnix系の環境（Linux、macOSなど）で使われることが多く、WindowsPathはWindows環境で使われることが多いです。

pathlibオブジェクトを文字列に変換するためには、str()関数を使用することができます。これにより、PosixPathやWindowsPathオブジェクトを必要な場所で文字列に変換することができます。

python
import pathlib

path = pathlib.Path("example.txt")
path_str = str(path)

一方で、pathlibはパスを直接操作するために/演算子をオーバーロードしています。これにより、パスオブジェクト同士を結合する際に直感的なコードを書くことができます。

python
import pathlib

path1 = pathlib.Path("dir1")
path2 = pathlib.Path("dir2")
combined_path = path1 / path2

したがって、pathlibライブラリはパスをより柔軟に操作できるために、環境に依存しないパスオブジェクトの型を返す設計となっています。必要に応じて、パスオブジェクトを文字列に変換することで、従来の文字列操作を行うこともできます。


WindowsやUNIXといったOSの違いに左右されないようにしているというのが回答のようですね。納得！！！
本題に戻ります。
A34-180316_GML.zipは解凍してみるとA34a-180316.shp(ポリゴンデータ)、A34b-180316.shp(ポイントデータ)、A34c-180316.shp(ラインデータ)の３つのシェープファイルを持っているようなので、それぞれべつのGeoDataFrameデータとして読み込みます。
import urllib.request
import zipfile

zip_url = "https://nlftp.mlit.go.jp/ksj/gml/data/A34/A34-17/A34-180316_GML.zip"

# ZIPファイルをダウンロード
urllib.request.urlretrieve(zip_url, "world_herirages.zip")

import pathlib
# ZIPファイルを解凍するフォルダー'data'を現在のディレクトリの下に作成
NOTEBOOK_PATH = pathlib.Path().resolve()
DATA_DIRECTORY = NOTEBOOK_PATH / "data"

# ZIPファイルを解凍して、'data'に保存
with zipfile.ZipFile("world_herirages.zip", "r") as zip_ref:
    zip_ref.extractall(DATA_DIRECTORY)

!pip install geopandas
import geopandas
shapefile_path = str(DATA_DIRECTORY) + "/A34-180316_GML/A34a-180316.shp"
world_heritages_poly = geopandas.read_file(shapefile_path)

shapefile_path = str(DATA_DIRECTORY) + "/A34-180316_GML/A34b-180316.shp"
world_heritages_point = geopandas.read_file(shapefile_path)

shapefile_path = str(DATA_DIRECTORY) + "/A34-180316_GML/A34c-180316.shp"
world_heritages_line = geopandas.read_file(shapefile_path)

world_heritages_poly.head()


A34a_001A34a_002A34a_003A34a_004A34a_005A34a_006A34a_007A34a_008A34a_009geometry
011–@—²Ž›’nˆæ‚Ì•§‹³Œš‘¢•¨]“Þ—ÇŒ§15.3–@—²Ž›ú@AúAAúCAúE199312POLYGON ((135.73529 34.61004, 135.73536 34.609...
121–@—²Ž›’nˆæ‚Ì•§‹³Œš‘¢•¨]“Þ—ÇŒ§15.3–@‹NŽ›ú@AúAAúCAúE199312POLYGON ((135.74658 34.62322, 135.74659 34.622...
232•P˜Hé]•ºŒÉŒ§107.0•P˜Héú@AúC199312POLYGON ((134.68559 34.83558, 134.68638 34.836...
343ŒÃ“s‹ž“s‚Ì•¶‰»à]Ž ‰êŒ§A‹ž“s•{1056.0‰ê–Εʗ‹_ŽÐúAAúC199412MULTIPOLYGON (((135.75183 35.06112, 135.75205 ...
453ŒÃ“s‹ž“s‚Ì•¶‰»à]Ž ‰êŒ§A‹ž“s•{1056.0‰ê–ÎŒä‘c_ŽÐúAAúC199412POLYGON ((135.77132 35.03817, 135.77165 35.039...

表示が文字化け！またまたChatGPT先生に質問。
User
world_heritages_point.head()　の表示が文字化けするのですが
ChatGPT

文字化けが発生する場合、GeoPandasのデータフレームで日本語が正しく表示されない可能性があります。文字エンコーディングを指定してデータを読み込むことで、文字化けを解決することができます。

GeoPandasでデータを読み込む際に、encodingパラメータを使用して文字エンコーディングを指定します。一般的に、日本語の文字エンコーディングには'utf-8'または'shift-jis'が使用されます。以下のようにencodingパラメータを指定してデータを読み込んでみてください。

python
world_heritages_point = geopandas.read_file(shapefile_path, encoding='shift-jis')
world_heritages_point.head()

上記の例では、shapefile_pathはデータのパスを指す変数であり、'shift-jis'という文字エンコーディングが使用されています。もし文字化けが解消されない場合、'utf-8'を試してみてください。

python
world_heritages_point = geopandas.read_file(shapefile_path, encoding='utf-8')
world_heritages_point.head()

適切な文字エンコーディングを指定することで、文字化けを解消できるはずです。ただし、データ自体に文字化けが含まれている場合は、それを解決するためには元のデータ自体を修正する必要があります。

ほうほう、ではまず’Shift-jis’でお試し。

A34a_001A34a_002A34a_003A34a_004A34a_005A34a_006A34a_007A34a_008A34a_009geometry
011法隆寺地域の仏教建造物‐奈良県15.3法隆寺@、A、C、E199312POLYGON ((135.73529 34.61004, 135.73536 34.609...
121法隆寺地域の仏教建造物‐奈良県15.3法起寺@、A、C、E199312POLYGON ((135.74658 34.62322, 135.74659 34.622...
232姫路城‐兵庫県107.0姫路城@、C199312POLYGON ((134.68559 34.83558, 134.68638 34.836...
343古都京都の文化財‐滋賀県、京都府1056.0賀茂別雷神社A、C199412MULTIPOLYGON (((135.75183 35.06112, 135.75205 ...
453古都京都の文化財‐滋賀県、京都府1056.0賀茂御祖神社A、C199412POLYGON ((135.77132 35.03817, 135.77165 35.039...

どうやら成功したようです。

これをとりあえずシンプルにplot()してみようと思います。
!pip install japanize_matplotlib
import matplotlib.pyplot as plt
import japanize_matplotlib

fig, ax = plt.subplots()

# ポリゴンデータのプロット
world_heritages_poly.plot(ax=ax, color='blue', label='地域')

# ポイントデータのプロット
world_heritages_point.plot(ax=ax, color='red', label='代表地点', markersize=2)

# ラインデータのプロット
world_heritages_line.plot(ax=ax, color='green', label='線')

# 凡例の作成
proxy_poly = plt.Rectangle((0, 0), 1, 1, fc='blue')
proxy_point = plt.Rectangle((0, 0), 1, 1, fc='red')
proxy_line = plt.Rectangle((0, 0), 1, 1, fc='green')

# 凡例の表示
ax.legend([proxy_poly, proxy_point, proxy_line], ['地域', '代表地点', '線'])

# プロットの表示
plt.show()

なんとなく雰囲気はわかりますが、日本地図を下敷きにしないと、なんのことやらさっぱりわからないですね。これをFoliumの地図上に表示したいと思います。
まずはfoliumの準備。適当な大きさの日本地図を表示します。
import folium

# Foliumの地図オブジェクトを作成
# 地図の中心は東京を指定
map = folium.Map(location=[35.6895, 139.6917], zoom_start=4)

map



インタラクティブな地図はblogerには簡単にアップロードできないので画像のコピペです。
この地図に先ほどの世界文化遺産に関するGeoDataFrameをプロットします。
import folium

# Foliumの地図オブジェクトを作成
# 地図の中心は東京を指定
map = folium.Map(location=[35.6895, 139.6917], zoom_start=4)

# ポリゴンをプロット
folium.GeoJson(world_heritages_poly).add_to(map)

# ポリゴンをプロット
folium.GeoJson(world_heritages_point).add_to(map)

# ラインストリングをプロット
folium.GeoJson(world_heritages_line).add_to(map)

map


これはpointにマーカーを、PolygonとLineは多角形や線（マーカーに隠れて見えません）を描画しています。アップするとそれらも見えます。わかりやすいように多角形(緑)と線(赤)に変更します。
import folium

# Foliumの地図オブジェクトを作成
# 地図の中心は東京を指定
map = folium.Map(location=[35.6895, 139.6917], zoom_start=4)

# ポリゴンをプロット
folium.GeoJson(
    world_heritages_poly,
    name='polygon',
    style_function=lambda feature: {
        'fillColor': 'green',
        'color': 'green',
        'weight': 2,
        'fillOpacity': 0.5
    }
).add_to(map)

# ポリゴンをプロット
folium.GeoJson(world_heritages_point).add_to(map)

# ラインストリングをプロット
folium.GeoJson(
    world_heritages_line,
    name='line',
    style_function=lambda x: {'color': 'red', 'weight': 2}
).add_to(map)
map



東海付近をアップしてみます。
京都・奈良周辺は当然多いとして、案外、富士山周辺にたくさんの世界文化遺産が存在していることに驚きました。でも、これではそれぞれの場所が何なのかがさっぱりわかりません。ぜひ、ポップアップを設定して、マーカーなどをクリックしたら、名称がでるようにしたいですね。
import folium

# Foliumの地図オブジェクトを作成
# 地図の中心は東京を指定
map = folium.Map(location=[35.6895, 139.6917], zoom_start=4)

# ポリゴンをプロット
folium.GeoJson(
    world_heritages_poly,
    name='polygon',
    style_function=lambda feature: {
        'fillColor': 'green',
        'color': 'green',
        'weight': 2,
        'fillOpacity': 0.5
    }
).add_to(map)

# ポイントをプロット
# GeoDataFrameから必要な列を取り出し、ポップアップに設定
for _, row in world_heritages_point.iterrows():
    popup_text = row['A34b_003']  # ポップアップに表示する列の値
    popup = folium.Popup(popup_text, max_width=200)
    folium.Marker([row.geometry.y, row.geometry.x], popup=popup).add_to(map)



# ラインストリングをプロット
folium.GeoJson(
    world_heritages_line,
    name='line',
    style_function=lambda x: {'color': 'red', 'weight': 2}
).add_to(map)
map



「富士山－信仰の対象と芸術の源泉」だそうです。多くのマーカーをクリックしましたが、全部同じポップアップでした。そうか、富士山か…。
しかし、四国に世界文化遺産がないというのは、とても寂しいですね。

これに四国88箇所霊場の位置をプロットして仕上げにしたいと思います。
#Google ColabでGoogle Driveを使用するときのおまじない
from google.colab import drive
drive.mount('/content/drive')

#四国巡礼88箇所の位置情報を記録したShapeファイルを読み取り
gdf_shikoku88 = gpd.read_file('/content/drive/MyDrive/Colab Notebooks/202303_Shikoku88Data/shikoku_88temples.shp')

for index, row in gdf_shikoku88.iterrows():
    popup_text = '第' + str(index) + '番霊場<br>' + row['Name']  # ポップアップに表示する列の値
    popup = folium.Popup(popup_text, max_width=200)
    folium.Marker(
        location=[row.geometry.y, row.geometry.x],
        icon=folium.Icon(icon='cloud', color='red'),  # マーカーの種類を変更
        popup=popup,
    ).add_to(map)

map

水色マーカーは世界文化遺産、赤色マーカーは四国88箇所霊場です！！！
本編はここで完結します。





















コメント

コメントを投稿