我收到了以下电子邮件,并希望确保每个人都能获得此问题的答案:
你好,
我想使用 cartopy 设置一个简单的纬度经度图,它穿过日期变更线,左侧显示东亚,右侧显示北美西部。以下谷歌地图大致是我所追求的:
https://maps.google.co.uk/?ll=56.559482,-175.253906&spn=47.333523,133.066406&t=m&z=4
这可以用 Cartopy 完成吗?
我收到了以下电子邮件,并希望确保每个人都能获得此问题的答案:
你好,
我想使用 cartopy 设置一个简单的纬度经度图,它穿过日期变更线,左侧显示东亚,右侧显示北美西部。以下谷歌地图大致是我所追求的:
https://maps.google.co.uk/?ll=56.559482,-175.253906&spn=47.333523,133.066406&t=m&z=4
这可以用 Cartopy 完成吗?
好问题。这可能会一次又一次地出现,所以在实际回答您的具体问题之前,我将逐步完成。为了将来参考,以下示例是使用 cartopy v0.5 编写的。
首先,重要的是要注意默认的“纬度经度”(或更技术上是 PlateCarree)投影在 -180 到 180 的正向范围内工作。这意味着您无法绘制超出此范围的标准 PlateCarree 投影。这有几个很好的原因,其中大部分归结为这样一个事实,即在投影矢量和栅格(例如简单的海岸线)时,cartopy 必须做更多的工作。不幸的是,您尝试制作的情节正是需要此功能。为了将此限制放入图片中,默认的 PlateCarree 投影如下所示:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
proj = ccrs.PlateCarree(central_longitude=0)
ax1 = plt.axes(projection=proj)
ax1.stock_img()
plt.title('Global')
plt.show()
您可以在此地图上绘制的任何单个矩形都可以合法地放大区域(这里有一些更高级的代码,但图片值 1000 字):
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
box = sgeom.box(minx=-90, maxx=45, miny=15, maxy=70)
x0, y0, x1, y1 = box.bounds
proj = ccrs.PlateCarree(central_longitude=0)
ax1 = plt.subplot(211, projection=proj)
ax1.stock_img()
ax1.add_geometries([box], proj, facecolor='coral',
edgecolor='black', alpha=0.5)
plt.title('Global')
ax2 = plt.subplot(212, projection=proj)
ax2.stock_img()
ax2.set_extent([x0, x1, y0, y1], proj)
plt.title('Zoomed in area')
plt.show()
不幸的是,您想要的绘图需要 2 个具有此投影的矩形:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
box = sgeom.box(minx=120, maxx=260, miny=15, maxy=80)
proj = ccrs.PlateCarree(central_longitude=0)
ax1 = plt.axes(projection=proj)
ax1.stock_img()
ax1.add_geometries([box], proj, facecolor='coral',
edgecolor='black', alpha=0.5)
plt.title('Target area')
plt.show()
因此,无法使用标准 PlateCarree 定义绘制跨越日期变更线的地图。相反,我们可以更改 PlateCarree 定义的中心经度,以允许在我们的目标区域中绘制一个框:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
box = sgeom.box(minx=120, maxx=260, miny=15, maxy=80)
x0, y0, x1, y1 = box.bounds
proj = ccrs.PlateCarree(central_longitude=180)
box_proj = ccrs.PlateCarree(central_longitude=0)
ax1 = plt.subplot(211, projection=proj)
ax1.stock_img()
ax1.add_geometries([box], box_proj, facecolor='coral',
edgecolor='black', alpha=0.5)
plt.title('Global')
ax2 = plt.subplot(212, projection=proj)
ax2.stock_img()
ax2.set_extent([x0, x1, y0, y1], box_proj)
plt.title('Zoomed in area')
plt.show()
希望这向您展示了实现目标地图必须做的事情,上面的代码可能有点复杂,以实现您的目标,所以稍微简化一下,我编写的用于生成您想要的绘图的代码将类似于:
import cartopy.feature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent([120, 260, 15, 80], crs=ccrs.PlateCarree())
# add some features to make the map a little more polished
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines('50m')
plt.show()
这是一个很长的答案,希望我不仅回答了这个问题,而且使地图制作和地图绘制的一些更复杂的细节更加清晰,以帮助解决您未来可能遇到的任何问题。
干杯,
有关上述benjimin评论的更多详细信息,
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.ticker import AutoMinorLocator, FixedLocator, MultipleLocator
def map_common(ax1,gl_loc=[True,True,False,True],gl_lon_info=range(-180,180,60),gl_dlat=30):
ax1.coastlines(color='silver',linewidth=1.)
gl = ax1.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=0.6, color='gray', alpha=0.5, linestyle='--')
gl.ylabels_left = gl_loc[0]
gl.ylabels_right = gl_loc[1]
gl.xlabels_top = gl_loc[2]
gl.xlabels_bottom = gl_loc[3]
gl.xlocator = FixedLocator(gl_lon_info)
gl.ylocator = MultipleLocator(gl_dlat)
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlabel_style = {'size': 11, 'color': 'k'}
gl.ylabel_style = {'size': 11, 'color': 'k'}
lon_boundary=np.arange(-240,-60,1.)
lat_boundary=np.arange(15,75,1.)
data=np.ones([lat_boundary.shape[0]-1,lon_boundary.shape[0]-1]) ## Data dimension is 1 less than boundaries
data=data*lat_boundary[:-1,None]
lon_offset=-150 ##
x,y=np.meshgrid(lon_boundary-lon_offset,lat_boundary)
fig=plt.figure()
fig.set_size_inches(7.5,5) ## (xsize, ysize)
ax1=fig.add_subplot(111,projection=ccrs.PlateCarree(central_longitude=lon_offset))
ax1.set_extent([-250,-50,10,80],crs=ccrs.PlateCarree())
props=dict(vmin=0,vmax=90,cmap=plt.cm.get_cmap('bone'),alpha=0.8)
cs=ax1.pcolormesh(x,y,data,**props)
ax1.set_title('Lon_Offset=-90')
map_common(ax1,gl_lon_info=[-180,-120,-60,120,],gl_dlat=15)
fnout='./map_over_dateline.png'
#plt.show()
plt.savefig(fnout,bbox_inches='tight',dpi=150)