やりたいこと
ラスタ画像上で任意地点の座標を指定して，値を拾うってのをやります．多少語弊があるかもしれませんが，QGISでいうところのPoint Sampling Toolみたいなやつですね．
今回は，国土数値情報の土地利用細分メッシュを用いて，指定した座標の土地利用の画素値をサンプリングしてみます．
使用ライブラリ
今回使うデータがGeoTiff形式のラスタデータなので，読み込みにみんな大好きGdalを使います．
インストールは以下のコマンドでできます．
pipの人
pip install GDAL
condaの人
conda install -c conda-forge gdal
使用データ
繰り返しになりますが，国土数値情報の土地利用細分メッシュを使います．今回は四国周辺の2014年度版のものを使ってます．
以下のリンクからダウンロードできます．
プログラム全体
from osgeo import gdal
import numpy as np
file_path = './input/landuse.tif'
ds = gdal.Open(file_path)
lu = ds.ReadAsArray()
Xsize = ds.RasterXSize
Ysize = ds.RasterYSize
x_min, x_res, rtn_x, y_max, rtn_y, y_res = ds.GetGeoTransform()
x_max = x_min + Xsize * x_res
y_min = y_max + Ysize * y_res
lon_tick = np.arange( x_min, x_max, x_res )
lat_tick = np.arange( y_max, y_min, y_res )
lon, lat = np.meshgrid(lon_tick, lat_tick)
p_lon = 133. + 32.9/60
p_lat = 33. + 34./60
dist = ( lon - p_lon )**2 + ( lat - p_lat )**2
mindist = np.min(dist)
p_index = np.where( dist == mindist )
p_lu = lu[p_index]
print(p_lu)
説明
値の拾い方
任意地点の座標値を起点に，全画素までの距離を測り，その大きさが最小になる地点の画素値を拾うというアイデアです．
ライブラリのインポート
from osgeo import gdal
import numpy as np
データの読み込み，必要な情報の取得
file_path = './input/landuse.tif'
ds = gdal.Open(file_path)
まず，gdalのOpenでGeoTiff画像を読み込みます．
lu = ds.ReadAsArray()
Xsize = ds.RasterXSize
Ysize = ds.RasterYSize
x_min, x_res, rtn_x, y_max, rtn_y, y_res = ds.GetGeoTransform()
上から順に，データの配列化，水平・鉛直方向の画素数，各種座標値を取得します．
GetGeoTransformでは，左から順に，始点端経度（最小経度），水平解像度，回転，始点端緯度（最大緯度），回転，鉛直解像度を取得できます．この時，鉛直解像度は負の値になることに注意してください．
x_max = x_min + Xsize * x_res
y_min = y_max + Ysize * y_res
lon_tick = np.arange( x_min, x_max, x_res )
lat_tick = np.arange( y_max, y_min, y_res )
lon, lat = np.meshgrid(lon_tick, lat_tick)
取得した座標値から，最大緯度，最小経度を求めます．
画像の４隅の座標がわかったので，これと縦横の解像度から，緯度，経度のメッシュを作成します．
値を拾う
座標の入力
p_lon = 133. + 32.9/60
p_lat = 33. + 34./60
任意地点から全画素までの距離を測る
dist = ( lon - p_lon )**2 + ( lat - p_lat )**2
中学数学でやる距離の公式です．ほんとは平方根取らないとちゃんとした距離にならないですが，大きさを比較するだけなので今回は平方根取ってません．
mindist = np.min(dist)
p_index = np.where( dist == mindist )
p_lu = lu[p_index]
print(p_lu)
>>> [70]
距離の最小値を求め，それに当たるインデックスを取得します．最後に，土地利用の画像上で取得したインデックスを指定することで，指定した座標の画素値を取得できます．
ちなみに今回は画素値として70という値が取得できました．対応表でチェックしてみると，「建物用地」みたいです．
注意点
指定した地点が隣り合う画素のど真ん中にある時など，距離の最小値が２つ以上ある時は困ります．特に今回のように画素値が離散的な値の時は，どの値を選べばいいか問題が発生することもあります．（なかなかそんなことは起きないけど）
ぼくの場合は画素値が連続値（気温とか標高値とか）の場合は，平均とったりしています．理屈上あり得る問題なので注意してください．