Python-地球科学-大气科学-可视化绘图系列(二)——利用basemap叠加地图,并添加白化效果(代码+示例)

时间:2020-01-14
本文章向大家介绍Python-地球科学-大气科学-可视化绘图系列(二)——利用basemap叠加地图,并添加白化效果(代码+示例),主要包括Python-地球科学-大气科学-可视化绘图系列(二)——利用basemap叠加地图,并添加白化效果(代码+示例)使用实例、应用技巧、基本知识点总结和需要注意事项,具有一定的参考价值,需要的朋友可以参考一下。

推荐参考链接:https://cloud.tencent.com/developer/article/1559267

 1 白化单图代码:
 2 import numpy as np
 3 import xarray as xr
 4 from mpl_toolkits.basemap import Basemap
 5 import matplotlib.pyplot as plt
 6 from matplotlib.patches import Polygon
 7 import matplotlib.patches as mpatches
 8 
 9 ds = xr.open_dataset('2019072300.006.nc')
10 t =  ds['value']
11 lons = ds.lon.data
12 lats = ds.lat.data
13 temp = xr.DataArray(t.data.T, coords=[lats,lons], dims=['latitude','longitude'])
14 # 创建画图空间
15 fig, ax = plt.subplots()
16 m = Basemap(projection='cyl',resolution='i',llcrnrlon=lons.min(),llcrnrlat=lats.min(),
17         urcrnrlon=lons.max(),urcrnrlat=lats.max(),lon_0=120.,lat_0=90)
18 Lon,Lat = np.meshgrid(lons[:],lats[:])
19 X,Y = m(Lon,Lat)
20 
21 shp_info3 = m.readshapefile("CHN_adm_shp\\CHN_adm3",'states',drawbounds=False,linewidth = 0.4,zorder=10)
22 for info, shp in zip(m.states_info, m.states):
23     proid = info['NAME_1']  # 可以用notepad打开CHN_adm1.csv文件,可以知道'NAME_1'代表各省的名称
24     if proid == 'Hubei':
25         poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.8)
26         ax.add_patch(poly)
27     else:
28         poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)
29         ax.add_patch(poly)
30 
31 #市一级底图
32 shp_info2 = m.readshapefile("CHN_adm_shp\\CHN_adm2",'states',drawbounds=False,linewidth = 0.4,zorder=10)
33 for info, shp in zip(m.states_info, m.states):
34     proid = info['NAME_1']
35     if proid == 'Hubei':
36         poly = Polygon(shp,facecolor='None',edgecolor='black', lw=0.8)
37         ax.add_patch(poly)
38 cs=m.contourf(Lon,Lat,t.data.T,np.arange(0,1.1,0.1),cmap='rainbow')
39 legend_elements = [mpatches.Patch(facecolor='#a50026',label='0-0.1'),
40                    mpatches.Patch(facecolor='#da362a',label='0.1-0.2'),
41                    mpatches.Patch(facecolor='#f67a49',label='0.2-0.3'),
42                    mpatches.Patch(facecolor='#fdbf6f',label='0.3-0.4'),
43                    mpatches.Patch(facecolor='#feeda1',label='0.4-0.5'),
44                    mpatches.Patch(facecolor='#ebf7a3',label='0.5-0.6'),
45                    mpatches.Patch(facecolor='#b7e075',label='0.6-0.7'),
46                    mpatches.Patch(facecolor='#75c465',label='0.7-0.8'),
47                    mpatches.Patch(facecolor='#249d53',label='0.8-0.9'),
48                    mpatches.Patch(facecolor='#006837',label='0.9-1')]
49 plt.legend(handles=legend_elements, loc='lower right',title='NDVI')
50 # cbar = plt.colorbar(cs,orientation='horizontal',label='Potential')

个例效果:

 1 多图代码:
 2 import numpy as np
 3 import xarray as xr
 4 from mpl_toolkits.basemap import Basemap
 5 import matplotlib.pyplot as plt
 6 from matplotlib.patches import Polygon
 7 import matplotlib.patches as mpatches
 8 
 9 ds = xr.open_dataset('2019072300.006.nc')
10 t =  ds['value']
11 lons = ds.lon.data
12 lats = ds.lat.data
13 temp = xr.DataArray(t.data.T, coords=[lats,lons], dims=['latitude','longitude'])
14 # 创建画图空间
15 fig, ax = plt.subplots(2,2,figsize=(15,15))
16 for ii in np.arange(4):
17     print('ii:::::',ii)
18     plt.sca(ax[np.unravel_index(ii,(2,2))]) # 获取一维索引在二维矩阵中的索引
19     m = Basemap(projection='cyl',resolution='i',llcrnrlon=lons.min(),llcrnrlat=lats.min(),
20         urcrnrlon=lons.max(),urcrnrlat=lats.max(),lon_0=120.,lat_0=90)
21     Lon,Lat = np.meshgrid(lons[:],lats[:])
22     X,Y = m(Lon,Lat)
23 
24     shp_info3 = m.readshapefile("CHN_adm_shp\\CHN_adm3",'states',drawbounds=False,linewidth = 0.4,zorder=10)
25     for info, shp in zip(m.states_info, m.states):
26         proid = info['NAME_1']  # 可以用notepad打开CHN_adm1.csv文件,可以知道'NAME_1'代表各省的名称
27         if proid == 'Hubei':
28             poly = Polygon(shp,facecolor='None',edgecolor='b', lw=0.8)
29             ax[np.unravel_index(ii,(2,2))].add_patch(poly)
30         else:
31             poly = Polygon(shp,facecolor='w',edgecolor='w', lw=0.8)
32             ax[np.unravel_index(ii,(2,2))].add_patch(poly)
33 
34     #市一级底图
35     shp_info2 = m.readshapefile("CHN_adm_shp\\CHN_adm2",'states',drawbounds=False,linewidth = 0.4,zorder=10)
36     for info, shp in zip(m.states_info, m.states):
37         proid = info['NAME_1']
38         if proid == 'Hubei':
39             poly = Polygon(shp,facecolor='None',edgecolor='black', lw=0.8)
40             ax[np.unravel_index(ii,(2,2))].add_patch(poly)
41 
42     cs=m.contourf(Lon,Lat,t.data.T,np.arange(0,1.1,0.1),cmap='rainbow')
43 
44 cbar = plt.colorbar(cs, ax=ax.ravel().tolist(), orientation='horizontal')
45 plt.savefig('multi_masked.jpg',bbox_inches = 'tight')

效果如下:

原文地址:https://www.cnblogs.com/zhanling/p/12193031.html