matplotlib のベースマップを使用してマップを描画しています。データは世界中に散らばっていますが、すべてのデータを大陸に残し、それらを海に落としたいだけです。データをフィルタリングする方法はありますか、またはデータをカバーするために海を再度描画する方法はありますか?
5 に答える
matplotlib.basemapには次のメソッドがあります。 is_land(xpt, ypt)
True
指定されたx、yポイント(投影座標内)が陸地上にある場合は戻ります。それ以外の場合は戻りますFalse
。土地の定義は、クラスインスタンスに関連付けられたGSHHS海岸線ポリゴンに基づいています。陸域内の湖の上のポイントは、陸地ポイントとしてカウントされません。
詳細については、こちらを参照してください。
is_land()
すべてのポリゴンをループして、土地かどうかを確認します。データサイズが大きい場合、非常に遅くなります。points_inside_poly()
from matplotlib を使用して、ポイントの配列をすばやく確認できます。これがコードです。lakepolygons
湖のポイントを削除したい場合は、自分で追加できます。
私のPCで100000ポイントを確認するのに2.7秒かかりました。速度を上げたい場合は、ポリゴンをビットマップに変換できますが、これを行うのは少し困難です。次のコードがデータセットに対して十分に高速でない場合は教えてください。
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.nxutils as nx
def points_in_polys(points, polys):
result = []
for poly in polys:
mask = nx.points_inside_poly(points, poly)
result.extend(points[mask])
points = points[~mask]
return np.array(result)
points = np.random.randint(0, 90, size=(100000, 2))
m = Basemap(projection='moll',lon_0=0,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
x, y = m(points[:,0], points[:,1])
loc = np.c_[x, y]
polys = [p.boundary for p in m.landpolygons]
land_loc = points_in_polys(loc, polys)
m.plot(land_loc[:, 0], land_loc[:, 1],'ro')
plt.show()
HYRY の回答は、matplotlib の新しいバージョンでは機能しません (nxutils は非推奨です)。動作する新しいバージョンを作成しました:
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.path import Path
import numpy as np
map = Basemap(projection='cyl', resolution='c')
lons = [0., 0., 16., 76.]
lats = [0., 41., 19., 51.]
x, y = map(lons, lats)
locations = np.c_[x, y]
polygons = [Path(p.boundary) for p in map.landpolygons]
result = np.zeros(len(locations), dtype=bool)
for polygon in polygons:
result += np.array(polygon.contains_points(locations))
print result