geopandas使ってみたメモ

Share Button

順々に増やしていく。Jupyterの中で練習のために作ったゴミみたいなデモファイルをそのまま消すのも憚られるので順に転記していくというだけ。

とりあえずArcGISのシェイプファイルを読んで図示したりする。
https://www.esrij.com/products/japan-shp/

jpn = gpd.read_file('japan_ver81.shp')
jpn.plot()


geopandasは結局のところpandasを継承しているので、操作感もpandasに慣れていれば簡単。
たとえば愛知県だけ引っ張ってくるならKENカラムが愛知県に一致するものだけ抽出したデータフレームに対してプロットをかければいい。

jpn[jpn['KEN'] == '愛知県'].plot()


さらに複数の条件で引っ張るなら、たとえば九州全域

area = ['福岡県', '宮崎県', '佐賀県', '熊本県', '鹿児島県', '長崎県', '大分県']
kyushu = jpn[jpn['KEN'].isin(area)]
kyushu.plot()


あとは都市のJISコードで指定したりもできる。墨田区のJISコードは13107なので…

jpn[jpn.JCODE == '13107'].plot()


地図をプロットした上にポイントを重ねるなら、

from geopandas.geoseries import Point

jpn = gpd.read_file('japan_ver81.shp')
tokyo23d = jpn[jpn.KEN.str.startswith('東京') & jpn.SIKUCHOSON.str.endswith('区')]
tokyo23d = tokyo23d.plot(color='white', edgecolor='black')

p = gpd.GeoSeries([Point(139.774722222, 35.682222222)])
p.plot(ax=tokyo23d, marker='o', color='red', markersize=5)

Point()関数に緯度経度を放り込んだものをGeoSeries()に配列で与えれあげればいい。


ある地点(今回も日本橋)がどこ(都内23区レベル)に属するのかを調べる(ジオコーディング)には,上のtokyo23dとpを使って

from geopandas.geoseries import Point

jpn = gpd.read_file('japan_ver81.shp')
tokyo23d = jpn[jpn.KEN.str.startswith('東京') & jpn.SIKUCHOSON.str.endswith('区')]

p = gpd.GeoDataFrame({
    'name': ['日本橋駅', '東京スカイツリー'], 
    'geometry': [Point(139.774722222, 35.682222222), Point(139.81083333333333, 35.71)]
})
p.crs = tokyo23d.crs

p_in_area = gpd.sjoin(p, tokyo23d, how="inner", op='intersects')
display(p_in_area)

結果を見ると…


ちなみにここまで直接的にGeoDataFrameを作ってましたが、通常のpandas.DataFrameからの変換もとりあえず仮データを作ってやってみる。longitudeとlatitudeが別のカラムとして格納されていることを想定して、
from shapely.geometry import Point

df = pd.DataFrame({
    'name': ['日本橋駅', '東京スカイツリー'], 
    'latitude': [35.682222222, 35.71],
    'longitude': [139.774722222, 139.81083333333333]
})
df['geometry'] = list(zip(df.longitude, df.latitude))
df['geometry'] = df.geometry.apply(Point)
gdf = gpd.GeoDataFrame(df, geometry='geometry')
手順として、GeoDataFrameにおける位置情報(geometryカラム)は基本的にlongitude, latitudeがshapelyのPoint関数でまとめられている。もちろんこれをPoint(long, lat)で放り込めばいいだけなので,longitude, latitudeの2つのカラムをタプルでひとつのカラムにまとめた上で、applyでPoint関数に通して,最後にGeoDataFrame()関数で変換しているだけ.その時にはどれがgeometryカラムなのか指定するのを忘れずに。

簡単ですね。時代も時代なのでもっと流行ってほしい。SASで触るのめんどい。


附録: 世田谷区の駅とその位置情報
DMSでしかデータ持ってこれなさそうだったから関数書いたんだけどよくみたら普通に10進数表記で取れたので途中でやめた。

def dms_to_Point(lat, long):
    latitude = lat[0] + lat[1] / 60 + lat[2] / 3600
    longitude = long[0] + long[1] / 60 + long[2] / 3600
    return Point(longitude, latitude)
sta = gpd.GeoDataFrame(pd.DataFrame({
    'name': [
        '二子玉川駅', '下北沢駅', '三軒茶屋', '成城学園前駅', '経堂駅',
        '駒沢大学駅', '千歳烏山駅', '用賀駅', '池尻大橋駅', '明大前駅',
        '千歳船橋駅', '桜新町駅', '等々力駅', '下高井戸駅', '桜上水駅',
        '喜多見駅',   '尾山台駅', '上野毛駅', '豪徳寺駅', '奥沢駅'
    ],
    'geometry': [
        dms_to_Point(lat=[35, 36, 42], long=[139, 37, 36]),
        dms_to_Point(lat=[35, 39, 41], long=[139, 40, 1]),
        dms_to_Point(lat=[35, 38, 37], long=[139, 40, 17]),
        dms_to_Point(lat=[35, 38, 24], long=[139, 35, 57]),
        dms_to_Point(lat=[35, 39, 4], long=[139, 38, 10]),
        Point(139.660971, 35.633048),
        Point(139.601015, 35.668181),
        Point(139.633933, 35.626469),
        Point(139.684580, 35.650629),
        Point(139.650554, 35.668349),
        Point(139.624640, 35.647656),
        Point(139.645025, 35.631874),
        Point(139.647946, 35.608355),
        Point(139.640939, 35.666136),
        Point(139.631751, 35.667622),
        Point(139.587422, 35.636684),
        Point(139.653890, 35.606967),
        Point(139.638844, 35.611994),
        Point(139.647361, 35.653739),
        Point(139.672264, 35.603988)
    ]
}), geometry='geometry')

style="display:block; text-align:center;"
data-ad-layout="in-article"
data-ad-format="fluid"
data-ad-client="ca-pub-3546003055292762"
data-ad-slot="5749192034">

Share Button

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

This site uses Akismet to reduce spam. Learn how your comment data is processed.