MeteoInfoLab脚本示例：一页多图

MeteoInfo

在MeteoInfoLab中实现一页多图可以用subplot（http://www.meteothink.org/docs ... s/plot/subplot.html）函数，但是如果一个subplot中有两个或者以上的坐标系（axes或axesm）时就不适合用subplot了，可以用axes（或者axesm）函数的position参数来控制子图的位置。这里演示制作一个4个不同物种（BC, PM2.5, SO2, NOx）的排放源清单图，图中有(a), (b), (c), (d)四个子图，每个子图还要有南海脚图，总共有8个坐标系（axesm）。

脚本程序：

1. spnames = ['BC','PM2.5','SO2','NOx']
   
2. nns = ['(a)','(b)','(c)','(d)']
   
3. basedir = 'D:/KeyData/PMMUL/analysis'
   
4. mapdir = 'D:/Temp/Map'
   
5. sectors = ['TRANSPORT','AIR','ENERGY','INDUSTRY','RESIDENTIAL','SHIPS']
   
6. 
   
7. #Figure
   
8. figure(figsize=[822, 421], newfig=False)
   
9. lworld = shaperead(os.path.join(mapdir, 'country1.shp'))
   
10. lchina = shaperead(os.path.join(mapdir, 'bou2_4p.shp'))
    
11. lsts = shaperead(os.path.join(basedir, 'shape/stations-1.shp'))
    
12. sc_layer = shaperead(os.path.join(mapdir, 'bou1_4l.shp'))
    
13. 
    
14. #Loop
    
15. pn = 1
    
16. col = 0
    
17. row = 0
    
18. for spname,nn in zip(spnames, nns):
    
19.     varname = 'emi_' + spname.lower()
    
20.     print varname
    
21.     datadir = 'D:/KeyData/Emission/EDGAR_HTAP/2010/' + spname
    
22.     i = 0
    
23.     for sector in sectors:
    
24.         f = addfile(os.path.join(datadir, 'edgar_HTAP_' + spname + '_emi_' + sector + '_2010.0.1x0.1.nc'))
    
25.         bc = f[varname][[15,55],[72,137]]
    
26.         if i == 0:
    
27.             tbc = bc
    
28.         else:
    
29.             tbc = tbc + bc
    
30.         i += 1
    
31.     # Plot
    
32.     if col == 2:
    
33.         row += 1
    
34.         col = 0
    
35.     print row
    
36.     print col
    
37.     axesm(position=[0.5*col+0.01,0.5*(1-row)+0.05,0.4,0.44], tickfontname='Tahoma', tickfontsize=10)
    
38.     geoshow(lworld)
    
39.     geoshow(lchina)
    
40.     geoshow(lsts, facecolor='b', size=2, labelfield='Stations', fontname='Tahoma', fontsize=10, yoffset=12)
    
41.     levs = arange(5.0E-12, 2.0E-10, 5.0E-12)
    
42.     if spname == 'BC':
    
43.         levs = arange(1.0E-12, 4.0E-11, 1.0E-12)
    
44.     cols = makecolors(len(levs) + 1, cmap='MPL_OrRd', reverse=False)
    
45.     cols[0] = 'w'
    
46.     layer = imshowm(tbc, levs, colors=cols)
    
47.     colorbar(layer, fontname='Tahoma', fontsize=10)
    
48.     xlim(72,137)
    
49.     ylim(15,55)
    
50.     text(75, 51, nn, fontname='Arial', fontsize=12, bold=True)
    
51.     yticks(arange(20, 60, 10))
    
52.     #Add south China Sea   
    
53.     axesm(position=[0.5*col+0.345-0.016*col,0.5*(1-row)+0.06,0.06,0.12], axison=False)
    
54.     geoshow(sc_layer)
    
55.     xlim(106, 123)
    
56.     ylim(2, 23)
    
57.     pn += 1
    
58.     col += 1
    
59. #savefig(os.path.join(basedir, 'Figure/emis_2010.pdf'))

spnames = ['BC','PM2.5','SO2','NOx']<br /> nns = ['(a)','(b)','(c)','(d)']<br /> basedir = 'D:/KeyData/PMMUL/analysis'<br /> mapdir = 'D:/Temp/Map'<br /> sectors = ['TRANSPORT','AIR','ENERGY','INDUSTRY','RESIDENTIAL','SHIPS']<br /> <br /> #Figure<br /> figure(figsize=[822, 421], newfig=False)<br /> lworld = shaperead(os.path.join(mapdir, 'country1.shp'))<br /> lchina = shaperead(os.path.join(mapdir, 'bou2_4p.shp'))<br /> lsts = shaperead(os.path.join(basedir, 'shape/stations-1.shp'))<br /> sc_layer = shaperead(os.path.join(mapdir, 'bou1_4l.shp'))<br /> <br /> #Loop<br /> pn = 1<br /> col = 0<br /> row = 0<br /> for spname,nn in zip(spnames, nns):<br />     varname = 'emi_' + spname.lower()<br />     print varname<br />     datadir = 'D:/KeyData/Emission/EDGAR_HTAP/2010/' + spname<br />     i = 0<br />     for sector in sectors:<br />         f = addfile(os.path.join(datadir, 'edgar_HTAP_' + spname + '_emi_' + sector + '_2010.0.1x0.1.nc'))<br />         bc = f[varname][[15,55],[72,137]]<br />         if i == 0:<br />             tbc = bc<br />         else:<br />             tbc = tbc + bc<br />         i += 1<br />     # Plot<br />     if col == 2:<br />         row += 1<br />         col = 0<br />     print row<br />     print col<br />     axesm(position=[0.5*col+0.01,0.5*(1-row)+0.05,0.4,0.44], tickfontname='Tahoma', tickfontsize=10)<br />     geoshow(lworld)<br />     geoshow(lchina)<br />     geoshow(lsts, facecolor='b', size=2, labelfield='Stations', fontname='Tahoma', fontsize=10, yoffset=12)<br />     levs = arange(5.0E-12, 2.0E-10, 5.0E-12)<br />     if spname == 'BC':<br />         levs = arange(1.0E-12, 4.0E-11, 1.0E-12)<br />     cols = makecolors(len(levs) + 1, cmap='MPL_OrRd', reverse=False)<br />     cols[0] = 'w'<br />     layer = imshowm(tbc, levs, colors=cols)<br />     colorbar(layer, fontname='Tahoma', fontsize=10)<br />     xlim(72,137)<br />     ylim(15,55)<br />     text(75, 51, nn, fontname='Arial', fontsize=12, bold=True)<br />     yticks(arange(20, 60, 10))<br />     #Add south China Sea    <br />     axesm(position=[0.5*col+0.345-0.016*col,0.5*(1-row)+0.06,0.06,0.12], axison=False)<br />     geoshow(sc_layer)<br />     xlim(106, 123)<br />     ylim(2, 23)<br />     pn += 1<br />     col += 1<br /> #savefig(os.path.join(basedir, 'Figure/emis_2010.pdf'))

aas314

好图学习了

梁的丰

期待已久的一页多图啊....

rjys1219

厉害，好看{:eb502:}

virl

厉害

wuxing3838

好图 学习了!

liuxiunuli

好图 学习了!，谢谢楼主

小哥8号

请问楼主怎么让南海填色和大陆一致

CHL1998

不好意思老师，打扰您一下，我在使用subplot函数进行一页多图的绘制时遇到了这种问题，能够整体分框，但是最后的数据图还是在图片中间。
麻烦您指导一下，我该怎么修改
另外，“axesm(newfig=False)”这个语句我也加上过，但是没有效果  

代码如下：

1. #read files
   
2. f1 = addfile('D:/MODIS/MOD08_M3/MOD08_M3.A2010152.061.2017313162026.hdf')
   
3. v_aod1 = f1['Aerosol_Optical_Depth_Land_Ocean_Mean_Mean']
   
4. v_aod1.setdimrev(0, True)
   
5. 
   
6. f2 = addfile('D:/MODIS/MOD08_M3/MOD08_M3.A2011152.061.2017325012251.hdf')
   
7. v_aod2 = f2['Aerosol_Optical_Depth_Land_Ocean_Mean_Mean']
   
8. v_aod2.setdimrev(0, True)
   
9. 
   
10. #get data
    
11. aod1 = v_aod1['30:45' , '105:125']
    
12. aod1[aod1 == -9.999000474927016] = nan
    
13. 
    
14. aod2 = v_aod2['30:45' , '105:125']
    
15. aod2[aod2 == -9.999000474927016] = nan
    
16. 
    
17. #plot
    
18. subplot(2,2,1)
    
19. levels = arange(0, 1.0, 0.02)
    
20. geoshow('cn_province', linecolor='gray')
    
21. geoshow('country', linecolor='gray')
    
22. layer = contourfm(aod1,levels,interpolate=True,cmap='NCV_jet')
    
23. colorbar(layer)
    
24. 
    
25. subplot(2,1,2)
    
26. levels = arange(0, 1.0, 0.02)
    
27. geoshow('cn_province', linecolor='gray')
    
28. geoshow('country', linecolor='gray')
    
29. layer = contourfm(aod1,levels,interpolate=True,cmap='NCV_jet')
    
30. colorbar(layer)

复制代码

MeteoInfo

CHL1998 发表于 2020-3-11 21:11
不好意思老师，打扰您一下，我在使用subplot函数进行一页多图的绘制时遇到了这种问题，能够整体分框，但是 ...

subplot 函数前两个数字不能变，subplot 函数中可以用 axestype 参数，不再需要 axesm 函数语句。

subplot(2, 2, 1, axestype='map')

subplot(2, 2, 2, axestype='map')