Python读取nc数据和绘图

zhengyunxia

学习学习，多谢分享

zhengyunxia

如果嫌风场箭头过于密集了要怎么设置？

呼啦啦

小白一个，已经下载pyproj和basemap两个包了，读取shp文件报错，是不是还需要安装什么
  File "C:\Users\qxt-4\Anaconda3\lib\site-packages\mpl_toolkits\basemap\__init__.py", line 2134, in readshapefile
    raise IOError('cannot locate %s.shp'%shapefile)
OSError: cannot locate D:\Program Files\python_workspace\ncep_download\.shp

平流层的萝卜

呼啦啦 发表于 2020-3-28 09:47
小白一个，已经下载pyproj和basemap两个包了，读取shp文件报错，是不是还需要安装什么
  File "C:%users\q ...

shp的文件如果是中文，比较麻烦。建议改成英文，删除奇奇怪怪的符号。

呼啦啦

平流层的萝卜 发表于 2020-3-30 08:36
shp的文件如果是中文，比较麻烦。建议改成英文，删除奇奇怪怪的符号。

您好，我看别的帖子，您有回复x轴为时间纵轴为高度的某一经纬度的剖面图，我想问下怎么画，手里的是FML的再分析资料，三维数据，是否是先转换为四维，然后在画，还是循环呢，有没有代码可以分享一下？

平流层的萝卜

呼啦啦 发表于 2020-4-4 17:36
您好，我看别的帖子，您有回复x轴为时间纵轴为高度的某一经纬度的剖面图，我想问下怎么画，手里的是FML的 ...

您好，其实原理是完全一样的，只不过“欺骗”电脑，让纵坐标Y轴从纬度变成高度。画的时候仍然是二维数据，一个x的坐标和一个Y的坐标，对应一个值。这样就形成了一个二维的场。
我之前写了一个比较长的脚本，部分功能是实现画剖面，贴在这里了：

1. # -*- coding: utf-8 -*-
   
2. import numpy as np
   
3. import netCDF4 as nc
   
4. # from mpl_toolkits.basemap import Basemap
   
5. import matplotlib.pyplot as plt
   
6. import time2page as tp
   
7. import thetase
   
8. import matplotlib as mpl
   
9. 
   
10. # print "fadsfasdfasd"
    
11. # exit()
    
12. # filesample=open(r'D:\liaoningrenying\dongbeilengwo\thin_reso_data\data\filenamelist.txt')
    
13. # filesample=np.loadtxt(r'h:\liaoningrenying\dongbeilengwo\thin_reso_data\data\filenamelist.txt',dtype='string',skiprows=15)
    
14. filesample=np.loadtxt(r'h:\liaoningrenying\dongbeilengwo\thin_reso_data\data\filenamelist_filtered_by_rain_trend_newest.txt',dtype='string',skiprows=4)
    
15. ######################下边是filenamelist.txt的内容，复制到此
    
16. '''
    
17. #注释：第一列为以时间命名的文件名，第二列为这个时刻在这个文件里是第几个时次。默认这个时间是降水开始的时间，但实际上需要更改，可另附一列。
    
18. #第三列是6h降水最强、范围最大的时刻，基本都是位于中部；第四列的意义：第二列+第四列的数，即为第三列时间在该nc文件里对应的时次；第五列是6h降水开始在辽宁出现的时刻，紧接着是降水出现在辽宁的哪里
    
19. #第七列是降水开始的时间相对于降水最长的时间，向前了几个时次
    
20. #06z08jun2011.nc 10 由于没有对应的降水数据，改为00z14may2012.nc 9 12051408
    
21. 00z10sep2008.nc 9 08091014 1 08090920 西部 -3
    
22. 00z14may2012.nc 9 12051408 0 12051320 西部 -2
    
23. 06z19apr2005.nc 10 05041920 1 05041902 西部 -3
    
24. 12z04nov2006.nc 11 06110420 0 06110414 西部 -1
    
25. 12z07jun2006.nc 11 06060714 -1 06060708 西部 -1
    
26. 18z03may2008.nc 12 08050408 1 08050320 西部 -2
    
27. 18z18may2006.nc 12 06051714 -6 06051708 西部 -1
    
28. '''
    
29. # print filesample
    
30. # exit()
    
31. dirthindata=r'h:\liaoningrenying\dongbeilengwo\thin_reso_data\data\\'
    
32. def pingjun(gaokong):
    
33.     return np.mean(gaokong,axis=0)
    
34. # lagt=-4
    
35. for lagt in [-4,-3,-2,-1,0,1,2,3,4]:
    
36.     # lagt=-1;
    
37.     print lagt
    
38.     if lagt==-4:timelag='-24h'
    
39.     if lagt==-3:timelag='-18h'
    
40.     if lagt==-2:timelag='-12h'
    
41.     if lagt==-1:timelag='-06h'
    
42.     if lagt==0:timelag='+00h'
    
43.     if lagt==1:timelag='+06h'
    
44.     if lagt==2:timelag='+12h'
    
45.     if lagt==3:timelag='+18h'
    
46.     if lagt==4:timelag='+24h'
    
47. 
    
48.     z,w,t,d,u,v,r,clwc,ciwc,cc=[],[],[],[],[],[],[],[],[],[]
    
49.     vor=[]
    
50.     for file_pagenum in filesample:
    
51.         filename,page_origin,heavyrain_time,delta_p,delta_p2=file_pagenum[0],file_pagenum[1],file_pagenum[2],file_pagenum[3],file_pagenum[6];
    
52.         page=int(page_origin)+int(delta_p)+int(delta_p2)    #这个page是6h降水开始出现的时刻
    
53.         # page=int(page_origin)+int(delta_p)      #这个page是6h降水最强的时刻
    
54.         page=page-1     ###表示20时发现有降水/降水达到最强，则14时开始降水/降水加强
    
55.         page=page-1     ####再减1是因为对应的第1页在list里对应第0个时间维，换算成list数
    
56.         page=int(page)+lagt
    
57.         # print page,page_origin,delta_p,lagt
    
58.         # print filename,page
    
59.         ncfile=nc.Dataset(dirthindata+filename,'r')
    
60.         z1,w1,t1,d1,u1,v1,r1,clwc1,ciwc1,cc1=ncfile['z'][page,:,:,:],ncfile['w'][page,:,:,:],ncfile['t'][page,:,:,:],ncfile['d'][page,:,:,:],ncfile['u'][page,:,:,:],\
    
61.                                             ncfile['v'][page,:,:,:],ncfile['r'][page,:,:,:],ncfile['clwc'][page,:,:,:],ncfile['ciwc'][page,:,:,:],ncfile['cc'][page,:,:,:]
    
62. 
    
63.         vor1=ncfile['vo'][page,:,:,:]
    
64.         # exit()
    
65.         lats,lons,levs=ncfile['latitude'][:],ncfile['longitude'][:],ncfile['level'][:]
    
66.         z.append(z1);w.append(w1);t.append(t1);d.append(d1);u.append(u1)
    
67.         v.append(v1);r.append(r1);clwc.append(clwc1);ciwc.append(ciwc1)
    
68.         cc.append(cc1);vor.append(vor1)
    
69.         # print ncfile
    
70.         ncfile.close()
    
71.     # print lons
    
72.     z,w,t,d,u,v,r,clwc,ciwc,cc,vor=map(pingjun,(z,w,t,d,u,v,r,clwc,ciwc,cc,vor)) ####z shape :23 ,201,201
    
73.     lon,lat=np.meshgrid(lons,lats)
    
74.     qvdivs,qxs,qys,qss=[],[],[],[]
    
75.     for i in range(0,23):
    
76.         qvdiv_i,qx_i,qy_i,qs_i=thetase.qvdiv((lon,lat,u[i,:,:],v[i,:,:],t[i,:,:],r[i,:,:],levs[i]))
    
77.         qvdivs.append(qvdiv_i);qxs.append(qx_i);qys.append(qy_i);qss.append(qs_i)
    
78.     qvdivs=np.array(qvdivs);qxs=np.array(qxs);qys=np.array(qys);qss=np.array(qss)
    
79.     # print np.array(qvdivs)*1e6
    
80.     # exit()
    
81. ######do slice cross-section#####
    
82.     # plt.figure(figsize=(1,1.2))
    
83.     fig,ax=plt.subplots(1,1)
    
84.     # plt.rcParams["figure.figsize"] = (10,14)
    
85.     # print "do cross-section:\n" \
    
86.     #       "First ,latitude=41.5,then longitude=",lon_cross_measure
    
87.     # nn=list(lats).index(41.5)
    
88.     # def latslice(array):
    
89.     #     a=array[:,nn,:]
    
90.     #     return a
    
91.     # z,w,d,u,v,r,clwc,ciwc=map(latslice,(z,w,d,u,v,r,clwc,ciwc))
    
92.     #
    
93.     # xx,zz=np.meshgrid(lons,levs)
    
94.     # contour=ax.contour(xx,zz,w,colors='k');contour.clabel(fmt='%6.4f',inline_spacing=0,fontsize=7)
    
95.     # ax.set_xlim(118,127)
    
96.     # plt.show()
    
97. 
    
98.     # lon_cross_measure=123
    
99.     # nn=list(lons).index(lon_cross_measure)
    
100.     # def lonslice(array):
     
101.     #     a=array[:,:,nn]
     
102.     #     return a
     
103.     # #
     
104.     # z,w,t,d,u,v,r,clwc,ciwc,cc,vor=map(lonslice,(z,w,t,d,u,v,r,clwc,ciwc,cc,vor))
     
105.     # qvdiv,qx,qy,qs=map(lonslice,(qvdivs,qxs,qys,qss))
     
106.     # yy,zz=np.meshgrid(lats,levs)
     
107. 
     
108.     print clwc
     
109.     # exit()
     
110.     lat_cross_measure=41
     
111.     nn=list(lats).index(lat_cross_measure)
     
112.     def latslice(array):
     
113.         a=array[:,nn,:]
     
114.         return a
     
115.     z,w,t,d,u,v,r,clwc,ciwc,cc,vor=map(latslice,(z,w,t,d,u,v,r,clwc,ciwc,cc,vor))
     
116.     qvdiv,qx,qy,qs=map(latslice,(qvdivs,qxs,qys,qss))
     
117. 
     
118.     yy,zz=np.meshgrid(lons,levs)
     
119.     print clwc
     
120.     # exit()
     
121. 
     
122.     ts,e,ei,ew=thetase.rh_t_p_to_thetase((r,t,zz)) #,norm=mpl.colors.Normalize(vmin=300, vmax=330)
     
123.     de=e-ei;ew_ei=ew-ei
     
124. 
     
125. 
     
126.     ##### 计算非绝热加热潜热
     
127.     Cp=1.005
     
128.     L=(2500-2.4*(t-273.15))*1e3
     
129.     dqsdp=np.gradient(qs,axis=0)/np.gradient(levs)[:,None]/100
     
130.     # print dqsdp
     
131.     # exit()
     
132.     theta=t*(1000/levs)[:,None]**0.285
     
133.     Hcs=-1.0/(Cp*t/theta)*L*w*dqsdp
     
134.     # np.savetxt('test_Q.txt',Hcs)
     
135.     # exit()
     
136.     ##### output binary of ts for grads verifying
     
137.     # np.float32((r[::-1,::-1],t[::-1,::-1],v[::-1,::-1],w[::-1,::-1])).tofile('thetase_output_by_python_for_grads_verification.grd')
     
138.     # # print type(r[1,1])
     
139.     # # print ts.shape
     
140.     # print yy,'\n',zz[::-1,1]
     
141.     # exit()
     
142.     # tshade=ax.contourf(yy,zz,ts,range(295,350,3),cmap=plt.get_cmap('jet'));plt.colorbar(tshade)
     
143.     # thetase_contour=ax.contour(yy,zz,ts,range(295,350,5),colors='gray');thetase_contour.clabel(fmt='%3d',inline_spacing=-10,fontsize=14)
     
144.     # rcontour=ax.contour(yy,zz,r,range(10,91,10),colors='blue');rcontour.clabel(fmt='%02d',inline_spacing=-5,fontsize=12)
     
145.     # deshade=ax.contourf(yy,zz,de,[0,0.1,0.2,0.3,0.4],cmap=plt.get_cmap('coolwarm'));plt.colorbar(deshade)
     
146.     decontour=ax.contour(yy,zz,de,[0,0.1,0.2,0.3,0.4],colors='k');decontour.clabel(inline_spacing=-10,fontsize=8)#,[0,0.1,0.2,0.3,0.4]
     
147.     # qvdivshade=ax.contourf(yy,zz,qvdiv*1e6,range(-50,21,10),cmap=plt.get_cmap('coolwarm'));colorbar_instance=plt.colorbar(qvdivshade)
     
148.     # print colorbar_instance.ColorbarBase
     
149.     # exit()
     
150.     # vor_contour=ax.contour(yy,zz,vor*1e5,range(-10,11,2),colors='brown');vor_contour.clabel(fmt='%3d')
     
151.     # Hcs_shade=ax.contourf(yy,zz,Hcs,range(-300,301,50),cmap='coolwarm',extend='both');#plt.colorbar(Hcs_shade,extendfrac=0.2,orientation='horizontal')
     
152.     # Hcs_contour=ax.contour(yy,zz,Hcs);Hcs_contour.clabel()
     
153.     clwc_shade=ax.contourf(yy,zz,clwc*1e5,range(3,16,3),cmap=plt.get_cmap('Blues'),extend='both');#plt.colorbar(clwc_shade,extendfrac=0.2,orientation='horizontal')       ####对于clwc
     
154.     # ciwc_shade=ax.contourf(yy,zz,ciwc*1e5,range(2,11,2),cmap=plt.get_cmap('Blues'),extend='both');plt.colorbar(ciwc_shade,extendfrac=0.2)    ####对于ciwc
     
155.     # clwc_contour=ax.contour(yy,zz,clwc*1e5,range(0,21,2),colors='green');clwc_contour.clabel(fmt='%4.0f',inline_spacing=-10,fontsize=12)
     
156.     # clwc_contour=ax.contour(yy,zz,ciwc*1e5,range(2,11,2),colors='k');clwc_contour.clabel(fmt='%4.0f',inline_spacing=-10,fontsize=12)
     
157.     # quiver=ax.quiver(yy[:,::4],zz[:,::4],(u[:,::4]),-w[:,::4]*100,color='k',headwidth=3,headlength=4,width=0.002)
     
158.     tscontour=ax.contour(yy,zz,t-273.15,[-15,-7,-4,0],colors='r');tscontour.clabel(fmt='%d$^\circ$C',inline_spacing=-5,fontsize=12)
     
159.     # tscontour=ax.contour(yy,zz,t,[273.15],colors='r');tscontour.clabel(fmt='%6.2f',inline_spacing=-10,fontsize=12)
     
160.     # decontour=ax.contour(yy,zz,de*100,range(0,51,5),colors='k');decontour.clabel(fmt='%2d' ,inline_spacing=-10,fontsize=8) #,(0,0.5,1.0,1.5)
     
161.     # decontour=ax.contour(yy,zz,de,[0],colors='gray');decontour.clabel(inline_spacing=-10,fontsize=8) #,(0,0.5,1.0,1.5) np.arange(0,0.1,0.01)
     
162.     vor[vor<0]=np.nan;v0contour=ax.contour(yy,zz,v*vor,[0],colors='brown',linewidths=4)
     
163. 
     
164.     # qvdivcontour=ax.contour(yy,zz,qvdiv*1e5,(-4,-3,-2,-1),colors='darkgreen',linewidths=1,linestyles='-');qvdivcontour.clabel(fmt='%3d' ,inline_spacing=-5,fontsize=8)
     
165.     # cc_contour=ax.contour(yy,zz,cc,np.arange(0.1,1.1,0.2),colors='k');cc_contour.clabel(fmt='%3.1f',inline_spacing=-5,fontsize=12)
     
166.     # wcontour=ax.contour(yy,zz,w,np.arange(-1.0,1.1,0.2),colors='brown');wcontour.clabel(fmt='%5.1f',inline_spacing=-10,fontsize=12)
     
167.     zcontour=ax.contour(yy,zz,z/9.8,(500,1000,2000,3000,4300,5500,6500),colors='orange',linestyles='--');zcontour.clabel(fmt='%4d',inline_spacing=0,fontsize=8)
     
168.     # ciwc_contour=ax.contour(yy,zz,ciwc*1e5,range(0,10),colors='purple');ciwc_contour.clabel(fmt='%4.0f',inline_spacing=-10,fontsize=12)
     
169.     # print np.shape(-w*100)
     
170.     # it=2
     
171.     # stream=ax.streamplot(yy,zz,v,-100.0*w,density=2)
     
172. 
     
173.     # stream=plt.streamplot(yy,zz,v,100*w,density=2)
     
174.      #ax.set_ylim(500,1000)
     
175.     # plt.xticks(range(30,56,5),[str(x)+u'°N' for x in range(30,56,5)],fontsize=13)
     
176.     plt.xticks(range(110,131,5),[str(x)+u'°E' for x in range(110,131,5)],fontsize=13);ax.set_xlim(110,130)
     
177.     plt.yticks((1000,925,900,850,800,700,600,500,400,300,200),fontsize=13)
     
178.     # plt.xlabel('Latitude',fontsize=15);plt.ylabel('Level',fontsize=15)
     
179.     # ax.set_xlim(38,44);#ax.set_ylim(200,1000)
     
180.     ax.set_ylim(300,1000)
     
181.     plt.gca().invert_yaxis()    #######把y坐标轴反转
     
182.     # plt.gca().invert_xaxis()
     
183.     plt.grid(True)
     
184. 
     
185.     # plt.title(u'Thetase cross-section plot along 123°E \n'+timelag,fontsize=18)
     
186.     # plt.savefig(r'H:\liaoningrenying\dongbeilengwo\\'+str(lat_cross_measure)+r'_cross_CLWC\\'+timelag+u'假相当位温_垂直环流_液水冰水含量_剖线'+str(lat_cross_measure)+'.png')
     
187.     plt.savefig(r'H:\liaoningrenying\dongbeilengwo\\'+str(lat_cross_measure)+r'_cross_CLWC\\'+timelag+u'假相当位温_垂直环流_液水冰水含量_剖线'+str(lat_cross_measure)+'.png')
     
188.     # plt.savefig(r'H:\liaoningrenying\dongbeilengwo\\'+str(lat_cross_measure)+r'_cross_CLWC\\'+timelag+u'假相当位温_垂直环流_CLWC_e-ei_剖线'+str(lat_cross_measure)+'.png')
     
189.     # plt.savefig(r'H:\liaoningrenying\dongbeilengwo\\'+str(lon_cross_measure)+r'_cross_CLWC\\'+timelag+u'假相当位温_垂直环流_液水冰水含量_剖线'+str(lon_cross_measure)+'.png',dpi=600)
     
190. 
     
191.     # print r'H:\liaoningrenying\dongbeilengwo\\'+str(lon_cross_measure)+r'_cross\\'+timelag+u'假相当位温_垂直环流_液水冰水含量_剖线'+str(lon_cross_measure)+'.png'
     
192. 
     
193.     # plt.savefig(r'H:\liaoningrenying\dongbeilengwo\\'+str(lat_cross_measure)+r'_cross_Q\\'+timelag+u'_Q_垂直环流_液水冰水含量_剖线'+str(lat_cross_measure)+'.png')
     
194. 
     
195.     # plt.show()
     
196.     plt.close()
     
197.     # break

复制代码

呼啦啦

平流层的萝卜 发表于 2020-4-6 17:56
您好，其实原理是完全一样的，只不过“欺骗”电脑，让纵坐标Y轴从纬度变成高度。画的时候仍然是二维数据 ...

好的，谢谢，我看看。

LeoMSII

感谢分享，学习了

hitwjk

赞一个吧 好样的

中间

我也要学，谢谢分享