【问题标题】:Fill oceans in basemap [duplicate]在底图中填充海洋[重复]
【发布时间】:2018-07-15 04:55:03
【问题描述】:

我正在尝试在matplotlib.Basemap 上绘制 1x1 度数据,并且我想用白色填充海洋。但是,为了让海洋的边界沿着matplotlib绘制的海岸线,白色海洋掩膜的分辨率应该比我数据的分辨率高很多。

找了半天,我尝试了两种可能的解决方案:

(1) maskoceans()is_land() 函数,但由于我的数据分辨率低于底图绘制的地图,因此边缘看起来不太好。我也不想将我的数据插值到更高分辨率。

(2) m.drawlsmask(),但由于无法分配 zorder,因此 pcolormesh 图总是覆盖蒙版。

这段代码

import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.basemap as bm

#Make data
lon = np.arange(0,360,1)
lat = np.arange(-90,91,1)
data = np.random.rand(len(lat),len(lon))

#Draw map
plt.figure()
m = bm.Basemap(resolution='i',projection='laea', width=1500000, height=2900000, lat_ts=60, lat_0=72, lon_0=319)
m.drawcoastlines(linewidth=1, color='white')
data, lon = bm.addcyclic(data,lon)
x,y = m(*np.meshgrid(lon,lat))
plt.pcolormesh(x,y,data)
plt.savefig('1.png',dpi=300)

生成此图像:

添加m.fillcontinents(color='white') 会生成以下图像,这是我需要的,但要填充海洋而不是陆地。

编辑

m.drawmapboundary(fill_color='lightblue') 也会填满土地,因此无法使用。

理想的结果是海洋是白色的,而我用plt.pcolormesh(x,y,data) 绘制的结果显示在陆地上。

【问题讨论】:

  • 我不认为你在这里问什么很清楚。您可以使用minimal reproducible example 使您的问题可重现,人们可以运行该minimal reproducible example,您可以使用不想要的结果的图片来解释问题。 m.drawmapboundary(fill_color='lightblue') 用浅蓝色填充海洋。目前尚不清楚为什么这在这里行不通。
  • 感谢您的建议,但m.drawmapboundary(fill_color='lightblue') 也填满了土地。
  • 那么期望的结果究竟是什么?请注意,底图中没有海洋,海洋只是没有被陆地覆盖的一切。也许您可以考虑到这一点来制定您的问题?
  • 理想的结果是海洋是白色的,而我用plt.pcolormesh(x,y,data) 绘制的图只显示在陆地上。所以与此相反:i.stack.imgur.com/w92pA.png,那里的土地是白色的,plt.pcolormesh(x,y,data) 绘制的内容只出现在海洋上。

标签: python matplotlib mask matplotlib-basemap


【解决方案1】:

我找到了一个更好的解决方案,它使用地图中海岸线定义的多边形来生成覆盖海洋区域的matplotlib.PathPatch。此解决方案具有更好的分辨率并且速度更快:

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

输出如下:

原解决方案

您可以在basemap.maskoceans 中使用分辨率更高的数组,以使分辨率适合大陆轮廓。之后,您可以反转掩码并将掩码数组绘制在数据之上。

不知何故,当我使用地图的全部范围(例如,经度从 -180 到 180 和纬度从 -90 到 90)时,我才让 basemap.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

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)


##producing a mask -- seems to only work with full coordinate limits
lons2 = np.linspace(-180,180,10000)
lats2 = np.linspace(-90,90,5000)
lon2, lat2 = np.meshgrid(lons2,lats2)
x2,y2 = m(lon2,lat2)
pseudo_data = np.ones_like(lon2)
masked = bm.maskoceans(lon2,lat2,pseudo_data)
masked.mask = ~masked.mask

##plotting the mask
cmap = colors.ListedColormap(['w'])
pcol = ax.pcolormesh(x2,y2,masked, cmap=cmap)

plt.show()

结果如下:

【讨论】:

  • 这和this question有什么不同?
  • @ImportanceOfBeingErnest 呃!可能不是。我怎么从来没有找到那些......?
  • 只是检查一下,但即使它不同,将解决方案添加到另一个问题可能更有意义,人们可以在一个地方找到所有选项。
  • @ImportanceOfBeingErnest 唯一的问题可能是数据的低分辨率,因此需要为掩码创建第二个更高分辨率的数组...
  • 可能与分辨率和完整地图的使用相关:gis.stackexchange.com/questions/236201/…
猜你喜欢
  • 1970-01-01
  • 2022-09-26
  • 2020-05-08
  • 1970-01-01
  • 2016-04-02
  • 2021-05-30
  • 2021-05-14
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多