求助 画低层数据时遮盖低于地形高度的部分

Noza220

如题，在画NCEP数据700hpa/850hpa水汽通量散度的时候想抹去小于地形高度的部分。
ncl中有mask函数，grades也听说过。因为MI已经写好代码了所以想用MI试试能不能做到。
我自己想的是最简单的那种，用位势高度减去地形高度，小于0的部分填色成白色。但是在实现的时候发现用contourm可以描绘出轮廓，但是用contourfm填色却不行。可能是我图层设置的问题..?
是不是只能用循环把要遮盖地区的数据写成缺省？或者用grades进行maskout数据？

我对python了解不多，希望各位老师和同学们多多指点！非常感谢！！
以下是代码，感谢王老师的水汽通量散度的程序.

1. #NCEP 6H 水汽通量散度
   
2. print 'Open data files...'
   
3. f_Air = addfile('E:/HELLO WORLD/DATA/NCEP/Water Vapor Flux/700hPa 9.3-5 air.nc')
   
4. f_Uwnd = addfile('E:/HELLO WORLD/DATA/NCEP/Water Vapor Flux/700hPa 9.3-5 uwind.nc')
   
5. f_Vwnd = addfile('E:/HELLO WORLD/DATA/NCEP/Water Vapor Flux/700hPa 9.3-5 vwind.nc')
   
6. f_Rhum = addfile('E:/HELLO WORLD/DATA/NCEP/Water Vapor Flux/700hPa 9.3-5 rh.nc')
   
7. 
   
8. f_hgt = addfile('E:/HELLO WORLD/DATA/NCEP/Water Vapor Flux/700hPa 9.3-5 hgt.nc') #700hpa位势高度
   
9. f_z = addfile('E:/HELLO WORLD/DATA/NCEP/hgt/NECP.hgt.sfc.nc') #地面位势高度
   
10. 
    
11. print 'Read data array...'
    
12. tidx = 8
    
13. t = f_Air.gettime(tidx)
    
14. Air = f_Air['air'][tidx,:,:,:]
    
15. Uwnd = f_Uwnd['uwnd'][tidx,:,:,:]
    
16. Vwnd = f_Vwnd['vwnd'][tidx,:,:,:]
    
17. Rhum = f_Rhum['rhum'][tidx,:,:,:]
    
18. hgt = f_hgt['hgt'][tidx,:,:,:]
    
19. z = f_z['hgt'][:,:,:]
    
20. 
    
21. # Calculate
    
22. print 'Calculate...'
    
23. speed = sqrt(Uwnd*Uwnd+Vwnd*Vwnd) #风场
    
24. 
    
25. prs = 700 #level
    
26. g = 9.8
    
27. es = 6.112*exp(17.67*(Air-273.16)/(Air-29.65)) #饱和水汽压
    
28. qs = 0.62197*es/(prs-0.378*es)*1000 #饱和比湿
    
29. q = qs*Rhum/100 #比湿
    
30. qhdivg = hdivg(q*Uwnd/g,q*Vwnd/g)*100000 #units = 10-5*g/(s•cm2•hPa)
    
31. 
    
32. blank = z-hgt
    
33. 
    
34. #Plot
    
35. print 'Plot...'
    
36. proj = projinfo(proj='lcc', lon_0=105, lat_1=25, lat_2=37)
    
37. axesm(projinfo=proj, axison=True, gridlabel=True, frameon=True)
    
38. mlayer = shaperead('E:/HELLO WORLD/MeteoInfo/map/country.shp')
    
39. lsichuan = shaperead('E:/HELLO WORLD/MeteoInfo/map/sichuan.shp')
    
40. geoshow(mlayer, edgecolor='black',size=2)
    
41. geoshow(lsichuan,edgecolor='black',size=2)
    
42. 
    
43. levs1 = [-3,-2.5,-2, -1.5, -1, -0.5,0, 0.5,1,1.5,2,2.5,3]
    
44. levs2 = [0,10000]
    
45. 
    
46. layer1 = contourfm(qhdivg,levs1, cmap='BlueRed') #水汽通量散度
    
47. layer2 = quiverm(Uwnd, Vwnd, speed,size=20,color='black') #风场
    
48. layer3 = contourm(blank,levs2,colors='white')#地形遮盖
    
49. 
    
50. title('Water Vapor Flux Divergency (' + t.strftime('%Y-%m-%d_%H:%M:%S') + ')',fontsize=20)
    
51. colorbar(layer1)
    
52. axism([81, 132, 9, 47])
    
53. #savefig('E:/HELLO WORLD/results/090500 700hPa WVP Divergency.png')
    

复制代码

这是出来的图，用layer3 = contoufm(blank,levs2,colors='white')的话就直接不显示了....不明白为什么啊，是因为不能用两个contoufm函数么？

                               
登录/注册后可看大图

问的问题可能非常小白，谢谢大家看帖。

MeteoInfo

参考下面的代码，只会将 blank 大于10000的等值区画成白色，小于10000的等值区的颜色是完全透明的（不会压盖下面图层的内容）。绘图函数的zorder参数也很重要，一个图中有两个contourfm语句生成的填色图层会互相压盖，blank图层要放在水汽通量散度图层之上。

levs2 = [10000]
cols = [(255,255,255,0), 'w']
layer3 = contourfm(blank, levs2, colors=cols, zorder=2)

Noza220

MeteoInfo 发表于 2019-5-22 23:08
参考下面的代码，只会将 blank 大于10000的等值区画成白色，小于10000的等值区的颜色是完全透明的（不会压 ...

因为想把blankd大于0的部分填白色，所以我改成了levs2 = [0]，然后zorder=2的话会压盖地图边界，所以改成=1。于是问题全部解决了。

谢谢王老师！！！！