僅在matplotlib的大陸上繪圖

Question

我正在使用matplotlib中的底圖繪製地圖。數據遍佈世界各地，但我只想保留大陸上的所有數據並將其放在海洋上。有沒有辦法可以過濾數據，還是有辦法再次繪製海洋來覆蓋數據？

Answer 1

有方法matplotlib.basemap：is_land(xpt, ypt)

它返回True如果給定的X，Y點（投影座標）是在陸地上，False否則。土地的定義基於與類實例相關的GSHHS海岸線多邊形。陸地地區內的湖泊點數不計爲陸地點。請參閱here。

Answer 2

is_land()將循環所有的多邊形來檢查它是否是土地。對於大數據量，它非常緩慢。您可以使用matplotlib中的points_inside_poly()快速檢查一組點。這是代碼。它不檢查lakepolygons，如果你想刪除湖泊中的點，你可以添加你的自我。

花了2.7秒在我的電腦上檢查100000個點。如果你想要更快的速度，你可以將多邊形轉換爲位圖，但這樣做有點困難。請告訴我，以下代碼對於您的數據集來說速度不夠快。

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() 

Answer 3

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 

Answer 4

最簡單的方法是使用底圖的maskoceans。

meshgrid和插值後：

from scipy.interpolate import griddata as gd 
from mpl_toolkits.basemap import Basemap, cm, maskoceans 
xi, yi = np.meshgrid(xi, yi) 
zi = gd((mlon, mlat), 
      scores, 
      (xi, yi), 
      method=grid_interpolation_method) 
#mask points on ocean 
data = maskoceans(xi, yi, zi) 
con = m.contourf(xi, yi, data, cmap=cm.GMT_red2green) 
#note instead of zi we have data now. 

Answer 5

我回答this question，當我被告知，這將是更好地張貼我的答案在這裏。基本上，我的解決方案提取用於繪製Basemap實例的海岸線的多邊形，並將這些多邊形與地圖的輪廓組合，以生成覆蓋地圖海洋區域的matplotlib.PathPatch。

如果數據比較粗糙並且不希望數據插值，這特別有用。在這種情況下，使用maskoceans會產生非常粗糙的海岸線輪廓，看起來不太好。

這裏是我張貼的答案其他問題同樣的例子：

from matplotlib import pyplot as plt 
from mpl_toolkits import basemap as bm 
from matplotlib import colors 
import numpy as np 
import numpy.ma as ma 
from matplotlib.patches import Path, PathPatch 

fig, ax = plt.subplots() 

lon_0 = 319 
lat_0 = 72 

##some fake data 
lons = np.linspace(lon_0-60,lon_0+60,10) 
lats = np.linspace(lat_0-15,lat_0+15,5) 
lon, lat = np.meshgrid(lons,lats) 
TOPO = np.sin(np.pi*lon/180)*np.exp(lat/90) 

m = bm.Basemap(resolution='i',projection='laea', width=1500000, height=2900000, lat_ts=60, lat_0=lat_0, lon_0=lon_0, ax = ax) 
m.drawcoastlines(linewidth=0.5) 

x,y = m(lon,lat) 
pcol = ax.pcolormesh(x,y,TOPO) 

##getting the limits of the map: 
x0,x1 = ax.get_xlim() 
y0,y1 = ax.get_ylim() 
map_edges = np.array([[x0,y0],[x1,y0],[x1,y1],[x0,y1]]) 

##getting all polygons used to draw the coastlines of the map 
polys = [p.boundary for p in m.landpolygons] 

##combining with map edges 
polys = [map_edges]+polys[:] 

##creating a PathPatch 
codes = [ 
    [Path.MOVETO] + [Path.LINETO for p in p[1:]] 
    for p in polys 
] 
polys_lin = [v for p in polys for v in p] 
codes_lin = [c for cs in codes for c in cs] 
path = Path(polys_lin, codes_lin) 
patch = PathPatch(path,facecolor='white', lw=0) 

##masking the data: 
ax.add_patch(patch) 

plt.show() 

這將產生以下情節：

希望這是有益的人:)