Yuta NakataのBlog

Python / AWS / ITについて役立つ情報を発信します

ベクトルデータ・ラスターデータから任意の緯度経度のデータを抽出する方法

課題背景

地理情報系のデータを扱う際、ベクトルデータとラスターデータがあります。 これらのデータは、度々巨大であることが多く、素直にデータを読み込むと、読み込み処理が重いという現象があります。

そこで、今回は、大きなデータであっても、データサイズの影響を受けることなく、スムーズにデータを抽出するコードを紹介します。

ラスターデータの場合

rasterioを使います。

import rasterio as rio

def get_mesh(infile: str, lat: float, lon: float, N=3):
    with rio.open(infile) as dataset:
        py, px = dataset.index(lon, lat)
        window = rio.windows.Window(px - N // 2, py - N // 2, N, N)
        dat = dataset.read(window=window)
    return dat

ベクトルデータの場合

geopandasを使う場合

import geopandas as gpd

lat, lon = 35, 135
point = Point(lon, lat)
df = gpd.read_file("sample.shp", bbox=(lon-1,  lat-1, lon+1, lat+1))
for index, polygon in enumerate(df['geometry']):
    if polygon.contains(point):
        return df.loc[index]

geopandasを使わない場合

geopandasは内部的に、gdalを使っているため、AWS Lambdaなどでは、依存関係が多く、Lambdaコンテナを使っていれる必要があります。

そこで、一度geopandasで読んだファイルをcsvで書き出してから、操作する方法があります。

import geopandas as gpd

gdf = gpd.read_file("sample.shp")

gdf.to_csv("result.csv", index=False)

このファイルをもとに、操作します。

import shapely

lat, lon = 35, 135

df = pd.read_csv('result.csv')
df['geometry'] = shapely.from_wkt(df['geometry'])

geom = df['geometry'].to_numpy()
point = shapely.points([lon], [lat])
tree = shapely.STRtree(geom)
ret = tree.query(point, 'within')
res = df.iloc[ret[1]]
print(res)