Cartopy绘制极射赤面投影

wxj96

最近需要画北半球的图，基于Python的Cartopy地图库，画北半球极射赤面投影。找了很多都没有相关的经验分享，准备自己写一个。


用到cartopy库的两个坐标函数，Orthographic和PlateCarree。Orthographic用于画底图，PlateCarree表示数据本身的投影方法。几乎所有的气象数据都是格点或者站点数据，基本都符合PlateCarree坐标。开始以为transform表示转化到某个投影，其实transform表示数据本身符合哪种投影，转化到地图投影上。


具体代码如下

1. import matplotlib.pyplot as plt
   
2. import matplotlib as mpl
   
3. import cartopy.crs as ccrs
   
4. import cartopy.feature as cfeature
   
5. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
   
6. from netCDF4 import Dataset
   
7. import numpy as np
   
8. 
   
9. def make_map(ax):
   
10.     projection = ccrs.Orthographic(central_longitude=-90,central_latitude=90)
    
11.     ax.set_global()
    
12.     ax.coastlines(linewidth=0.5)
    
13.     '''标注坐标轴'''
    
14.     # ax.set_xticks(np.linspace(-180, 180, 1), crs=projection)
    
15.     # ax.set_yticks(np.linspace(-90, 90,5), crs=projection)
    
16.     '''zero_direction_label=True 有度的标识，False则去掉'''
    
17.     lon_formatter = LongitudeFormatter(zero_direction_label=True)
    
18.     lat_formatter = LatitudeFormatter()
    
19.     # ax.xaxis.set_major_formatter(lon_formatter)
    
20.     # ax.yaxis.set_major_formatter(lat_formatter)
    
21.     '''添加网格线'''
    
22.     #gl = ax.gridlines()
    
23.     #ax.grid()
    
24.     return ax
    
25. 
    
26. 
    
27. data = Dataset('hgt500_correlation_mon_mean.nc')
    
28. hgt = data.variables['correlation'][:,:,:,:]
    
29. lat = data.variables['lat'][:]
    
30. lon = data.variables['lon'][:]
    
31. data_proj = ccrs.PlateCarree(central_longitude=0)
    
32. fig = plt.figure()
    
33. proj = ccrs.Orthographic(central_longitude=-90,central_latitude=90)
    
34. ax = fig.add_subplot(projection = proj)
    
35. ax = make_map(ax)
    
36. ax.gridlines(color='gray', linestyle='--', xlocs = np.arange(0,390,20), ylocs = np.linspace(-80,80,9))
    
37. cf = ax.contour(lon, lat, hgt[0,0,:,:], colors='k',transform=data_proj)
    
38. plt.clabel(cf,fmt = '%5.2f')
    
39. plt.show()
    
40. 
    

复制代码

我不是卢姥爷

{:5_213:}{:5_213:}{:5_213:}{:5_213:}{:5_213:}

我不是卢姥爷

真的太牛逼了！

易烊启迪

博主您好，如果不画整个北半球，只画纬度为40-90的，经度为0-360的地图可以怎么修改代码吗？

wxj96

易烊启迪 发表于 2022-8-22 17:15
博主您好，如果不画整个北半球，只画纬度为40-90的，经度为0-360的地图可以怎么修改代码吗？

xlim和ylim用于选择绘图范围

yyh61669641

很好，就是正面中心是北美洲，改成正面是110E或120E才是需要的