Python. 在basemap中绘制矩形
我需要在我的底图上添加几个矩形。我需要四个矩形,它们的经纬度范围如下:
1) 左下角经度=-10,右上角经度=10,左下角纬度=35,右上角纬度=60
2) 左下角经度=10.5,右上角经度=35,左下角纬度=35,右上角纬度=60
3) 左下角经度=35.5,右上角经度=52,左下角纬度=30,右上角纬度=55
4) 左下角经度=-20,右上角经度=35,左下角纬度=20,右上角纬度=34.5
我的脚本在下面。我找到了“polygon”包来添加线条,但我不太清楚该怎么做。请帮帮我!!非常感谢你的帮助!
from mpl_toolkits.basemap import Basemap
m=basemaputpart.Basemap(llcrnrlon=-60, llcrnrlat=20, urcrnrlon=60, urcrnrlat=70, resolution='i', projection='cyl', lon_0=0, lat_0=45)
lon1=np.array([[-180.+j*0.5 for j in range(721)] for i in range(181)])
lat1=np.array([[i*0.5 for j in range(721)] for i in range(181) ])
Nx1,Ny1=m(lon1,lat1,inverse=False)
toplot=data[:,:]
toplot[data==0]=np.nan
toplot=np.ma.masked_invalid(toplot)
plt.pcolor(Nx1,Ny1,np.log(toplot),vmin=0, vmax=5)
cbar=plt.colorbar()
m.drawcoastlines(zorder=2)
m.drawcountries(zorder=2)
llcrnrlon = -10
urcrnrlon = 10
llcrnrlat = 35
urcrnrlat = 60
lower_left = (llcrnrlon, llcrnrlat)
lower_right= (urcrnrlon, llcrnrlat)
upper_left = (llcrnrlon, urcrnrlat)
upper_right= (urcrnrlon, urcrnrlat)
plot_rec(m, lower_left, upper_left, lower_right, upper_right)
然后我看到“类型错误:'tuple'对象不可调用”
我在这部分中,先添加了你建议的内容。
..
m.drawcoastlines(zorder=2)
m.drawcountries(zorder=2)
def plot_rec(m, lower_left, upper_left, lower_right, upper_right):
xs = [lower_left[-10], upper_left[-10],
lower_right[10], upper_right[10]]
ys = [lower_left[35], upper_left[60],
lower_right[35], upper_right[60]]
m.plot(xs,ys,latlon=True)
plt.show()
然后我在我的图中看不到任何框。我是不是需要放另一个,而不是 plt.show()??
另外,你能告诉我怎么在框里放数字吗(例如,在框的左上角放1)?怎么计算我所有数据点的值的总和,并计算(框内值的总和)与(所有数据点的值的总和)之间的百分比?我问得太多了……只要告诉我你能提供的内容就好,反正都很棒!!!
非常感谢!!
3 个回答
如果你想在地图上找那些边界沿着经线和纬线的矩形区域,可以使用下面这个函数。这里的 m 是基础地图对象。
def draw_screen_poly( minlat, maxlat, minlon, maxlon, m):
lons=np.hstack((np.repeat(minlon,10),\
np.linspace(minlon,maxlon, num=10),\
np.repeat(maxlon,10),\
np.linspace(maxlon,minlon, num=10)))
lats=np.hstack((np.linspace(minlat,maxlat, num=10),\
np.repeat(maxlat,10),\
np.linspace(maxlat,minlat, num=10),
np.repeat(minlat,10)))
m.plot(y=lats,x=lons,latlon=True, lw=2, color='navy', alpha=0.8)
x, y = m(lons, lats)
xy = zip(x,y)
poly = Polygon( np.asarray(xy), linewidth=3)
ax.add_patch(poly)
下面这个函数应该能满足你的需求。bmap 是你的基础地图对象,而 lonmin, lonmax, latmin, latmax 则定义了你想要显示的区域的经纬度范围。在你用 Basemap(...) 生成了 bmap 之后,你需要调用 plot_rectangle
函数。
def plot_rectangle(bmap, lonmin,lonmax,latmin,latmax):
xs = [lonmin,lonmax,lonmax,lonmin,lonmin]
ys = [latmin,latmin,latmax,latmax,latmin]
bmap.plot(xs, ys,latlon = True)
这件事比较棘手,因为在很多地图投影类型中,“矩形”并不是真正的“矩形”。所以当你说“矩形”时,你是指地图上的实际矩形,还是仅仅是像素空间中的矩形?这两者的处理方式是完全不同的。
但我们假设你想要的是地图上的矩形。最快的方法就是使用Basemap的plot方法,像这样:
def plot_rec(bmap, lower_left, upper_left, lower_right, upper_right):
xs = [lower_left[0], upper_left[0],
lower_right[0], upper_right[0],
lower_left[0], lower_right[0],
upper_left[0], upper_right[0]]
ys = [lower_left[1], upper_left[1],
lower_right[1], upper_right[1],
lower_left[1], lower_right[1],
upper_left[1], upper_right[1]]
bmap.plot(xs, ys, latlon = True)
这里的bmap
是你的地图,而lower_left
等则是这些角落的经纬度坐标。
更新:使用示例
你问了一个使用示例,所以这里有一个:
m=basemaputpart.Basemap(llcrnrlon=-60, llcrnrlat=20, urcrnrlon=60, urcrnrlat=70, resolution='i', projection='cyl', lon_0=0, lat_0=45)
# your other setting up the map code here
# here I draw the first rectangle
llcrnrlon = -10
urcrnrlon = 10
llcrnrlat = 35
urcrnrlat = 60
lower_left = (llcrnrlon, llcrnrlat)
lower_right= (urcrnrlon, llcrnrlat)
upper_left = (llcrnrlon, urcrnrlat)
upper_right= (urcrnrlon, urcrnrlat)
plot_rec(m, lower_left, upper_left, lower_right, upper_right) # This calls the function I defined before
# Rinse and repeat for the other lat/lon combos
plt.show()
你当然可以用列表推导式更优雅地生成正确的角点集合,但这个示例应该能让你入门。
更新 2
看起来这里有些混淆。plot_rec
是一个函数。它应该放在你的脚本中某个不与其他代码行在一起的地方。单独放着它是没有任何作用的。只有在你调用它的时候,它才会起作用,如下所示:
...
upper_left = (llcrnrlon, urcrnrlat)
upper_right= (urcrnrlon, urcrnrlat)
# This next line is where we call the function
plot_rec(m, lower_left, upper_left, lower_right, upper_right)