【OLD】在地圖上畫crosssection

比如我想沿著黑直線這個位置畫crosssection

首先，先畫這個條在地圖上，畫出上面這樣的圖

def seearea():  ##--read data and analysis--##  filename="4draw_950_lljandnonllj_uvz"  da=nc.Dataset(filename)  u=da.variables["u"][:]  v=da.variables["v"][:]  z=da.variables["z"][:]  lon=da.variables["lon"][:]  lat=da.variables["lat"][:]  da.close()  da=nc.Dataset("HGT.nc")  hgt=da.variables["hgt"][:]  da.close()  np.set_printoptions(threshold="nan")  wind=np.sqrt(u**2+v**2)  print(np.shape(wind))  plt.figure()  lon1=100  lon2=122  lat1=6  lat2=25  ##--draw map--##  m = Basemap(projection='rotpole',o_lat_p=84,o_lon_p=00,lon_0=0,lat_0=0,lat_1=lat1,lat_2=lat2,lon_1=lon1+1,lon_2=lon2-1,llcrnrlat=lat1,urcrnrlat=lat2,llcrnrlon=lon1,urcrnrlon=lon2,resolution='c')  m.drawcoastlines(color='k',linewidth=1)  m.drawcountries(color='k',linewidth=1)#2e5244  m.drawcoastlines(color='r',linewidth=0.5)  m.drawcountries(color='r',linewidth=0.5)#2e5244  m.drawparallels(np.arange(0.,90.,5),labels=[1,0,0,0])  m.drawmeridians(np.arange(0.,180,5),labels=[0,0,0,1])  plt.title("u,v,z,in 950hPa",size='x-large')  x, y = m(lon,lat)  ##--contour--##  csslev = np.arange(math.floor(np.min(z)),math.ceil(np.max(z)),4)  css = m.contour(x,y,z,csslev,linewidths=1,colors='lime')  css = m.contour(x,y,z,csslev,linewidths=0.5,colors='k')  plt.clabel(css, inline=0,fontsize=8)  plt.figtext(0.5,0.94,'height contour= %d to %d(unit: m)'%(math.floor(np.min(z)),math.ceil(np.max(z))),fontsize=12)  ##--contourf--##  cslev = np.arange(0,13,1)  #cslev = np.arange(0,1200,12)  cs = m.contourf(x,y,wind,cslev,cmap=cm.coolwarm,extend='both')  cs.cmap.set_over('k')  cs.cmap.set_under('b')  cbar = m.colorbar(cs,location='bottom',pad="8%")  cbar.set_label('wind speed (m/s)')  ##--vector--##  vectorsize=5  rho=7  Q = m.quiver(x[::rho,::rho],y[::rho,::rho],u[::rho,::rho],v[::rho,::rho],width=vectorsize,headwidth=vectorsize*2/3,headlength=vectorsize/2,headaxislength=vectorsize/2,minlength=0,units='dots',angles='xy', scale_units='xy', scale=30,minshaft=0.5)#30  qk = plt.quiverkey(Q, 0.1,1.1, 20, '20 m/s', labelpos='S')  ##--draw plot--##  #lat001=range(100,132)  #lon001=np.floor(np.arange(180,228,1.5))  theta=0.2-np.pi/2#0.875-np.pi/2  '''  da=nc.Dataset("C:\\Users\\kong\\PycharmProjects\\southcnseapycode\\crosssection\\4draw_crosssection")  lat001=da.variables["lat001"][:]  lon001=da.variables["lon001"][:]  da.close()  for ii in range(0,np.shape(lat001)[0]):    m.plot(x[lon001[ii],lat001[ii]],y[lon001[ii],lat001[ii]],'ko',markersize=2)    ii=ii+1  '''  lon001=np.arange(80,170)  lat001=((lon001-125)/math.tan(theta)+130)  #lat001=np.arange(120,140,0.1)  #lon001=np.floor((lat001-130)*math.tan(theta)+125)  for ii in range(0,np.shape(lat001)[0]):    meow=lat001[ii]    meow1=np.floor(lat001[ii])    meow2=np.ceil(lat001[ii])    if meow1==meow2:      m.plot(x[lon001[ii],lat001[ii]],y[lon001[ii],lat001[ii]],'ko',markersize=1)    else:      m.plot((x[lon001[ii],meow2]-x[lon001[ii],meow1])/(meow2-meow1)*(meow-meow1)+x[lon001[ii],meow1],(y[lon001[ii],meow2]-y[lon001[ii],meow1])/(meow2-meow1)*(meow-meow1)+y[lon001[ii],meow1],'ko',markersize=1)    ii=ii+1  #m.plot(x[210,120],y[210,120],'ro',markersize=5)  #print(lon[210,120],lat[210,120])  #m.plot(x[125,130],y[125,130],'ro',markersize=5)  u1=5  csslev = np.arange(50,1600,100) #######  css = m.contour(x,y,hgt,csslev,linewidths=0.3,colors='k')#0.5  #plt.quiver(x[210,110],y[210,110],u1*np.cos(theta),u1*np.sin(theta), scale=30)  ##--save image--##  plt.savefig("C:\\Users\\kong\\PycharmProjects\\southcnseapycode\\crosssection\\seeareaforcrosssection.png", dpi=500)#500#max dpi2300  #plt.show()  print(lon001)  print(lat001)

算lon, lat值時用到線性插值，其實不用也可以，不過因為這個crosssection與經線的夾角太小了，因此不用插值而用floor或ceil等來取整畫出來的圖會很丑。

用上面的代碼找出你想畫的地方以後﹐可以用下面的代碼把crosssection的數據輸出

import numpy as npimport netCDF4 as ncimport mathprint("IMPORT FINISH")##--read data--####--u--##da=nc.Dataset("4draw_average_u_3d")u=da.variables["u"][:]lon=da.variables["lon"][:]lat=da.variables["lat"][:]da.close()##--v--##da=nc.Dataset("4draw_average_v_3d")v=da.variables["v"][:]da.close()##--w--##da=nc.Dataset("4draw_average_w_3d")w=da.variables["w"][:]da.close()##--z--##da=nc.Dataset("4draw_average_z_3d")z=da.variables["z"][:]da.close()##--llj--##da=nc.Dataset("4draw_average_llj_3d")llj=da.variables["llj"][:]da.close()print("READ FINISH")##--line--##theta=0.2-np.pi/2lon001=np.arange(80,170)lat001=((lon001-125)/math.tan(theta)+130)lon001=list(lon001)lat001=list(lat001)if np.shape(lat001)[0]==np.shape(lon001)[0]:  ##--output--##  output_lon=np.zeros((np.shape(lat001)[0]))  output_lat=np.zeros((np.shape(lat001)[0]))  output_u=np.zeros((25,np.shape(u)[1],np.shape(lat001)[0]))  output_v=np.zeros((25,np.shape(u)[1],np.shape(lat001)[0]))  output_w=np.zeros((25,np.shape(u)[1],np.shape(lat001)[0]))  output_z=np.zeros((25,np.shape(u)[1],np.shape(lat001)[0]))  output_llj=np.zeros((25,np.shape(u)[1],np.shape(lat001)[0]))  lon=np.array(lon)    lat=np.array(lat)  i_out=0  for i_lat in lat001:    lat1=(np.floor(i_lat))    lat2=(np.ceil(i_lat))    if i_lat==lat2:      output_lon[i_out]=lon[lon001[i_out],i_lat]      output_lat[i_out]=lat[lon001[i_out],i_lat]      output_u[:,:,i_out]=u[:,:,lon001[i_out],i_lat]      output_v[:,:,i_out]=v[:,:,lon001[i_out],i_lat]      output_w[:,:,i_out]=w[:,:,lon001[i_out],i_lat]      output_z[:,:,i_out]=z[:,:,lon001[i_out],i_lat]      output_llj[:,:,i_out]=llj[:,:,lon001[i_out],i_lat]    else:      output_lon[i_out]=(lon[lon001[i_out],lat2]-lon[lon001[i_out],lat1])/(lat2-lat1)*(i_lat-lat1)+lon[lon001[i_out],lat1]      output_lat[i_out]=(lat[lon001[i_out],lat2]-lat[lon001[i_out],lat1])/(lat2-lat1)*(i_lat-lat1)+lat[lon001[i_out],lat1]      output_u[:,:,i_out]=(u[:,:,lon001[i_out],lat2]-u[:,:,lon001[i_out],lat1])/(lat2-lat1)*(i_lat-lat1)+u[:,:,lon001[i_out],lat1]      output_v[:,:,i_out]=(v[:,:,lon001[i_out],lat2]-v[:,:,lon001[i_out],lat1])/(lat2-lat1)*(i_lat-lat1)+v[:,:,lon001[i_out],lat1]      output_w[:,:,i_out]=(w[:,:,lon001[i_out],lat2]-w[:,:,lon001[i_out],lat1])/(lat2-lat1)*(i_lat-lat1)+w[:,:,lon001[i_out],lat1]      output_z[:,:,i_out]=(z[:,:,lon001[i_out],lat2]-z[:,:,lon001[i_out],lat1])/(lat2-lat1)*(i_lat-lat1)+z[:,:,lon001[i_out],lat1]      output_llj[:,:,i_out]=(llj[:,:,lon001[i_out],lat2]-llj[:,:,lon001[i_out],lat1])/(lat2-lat1)*(i_lat-lat1)+llj[:,:,lon001[i_out],lat1]    i_out=i_out+1  ##--save data--##  da=nc.Dataset("4draw_crosssection_125_130","w",format="NETCDF4")  da.createDimension("heisize",np.shape(u)[1])  da.createDimension("xsize",np.shape(lat001)[0])  da.createDimension("timesize",np.shape(u)[0])  da.createDimension("onesize",1)  da.createVariable("u","f8",("timesize","heisize","xsize"))  da.createVariable("v","f8",("timesize","heisize","xsize"))  da.createVariable("w","f8",("timesize","heisize","xsize"))  da.createVariable("z","f8",("timesize","heisize","xsize"))  da.createVariable("llj","f8",("timesize","heisize","xsize"))  da.createVariable("lon","f8",("xsize",))  da.createVariable("lat","f8",("xsize",))  da.createVariable("lon001","f8",("xsize",))  da.createVariable("lat001","f8",("xsize",))  da.createVariable("theta","f8",("onesize",))  da.variables["u"][:]=output_u  da.variables["v"][:]=output_v  da.variables["w"][:]=output_w  da.variables["z"][:]=output_z  da.variables["llj"][:]=output_llj  da.variables["lat"][:]=output_lat  da.variables["lon"][:]=output_lon  da.variables["lat001"][:]=lat001  da.variables["lon001"][:]=lon001  da.variables["theta"][:]=theta  da.description="crosssection"  da.author="konghoiio"  da.createdate="2016-11-08"  da.close()else:  print("error happened")print("ALL FINISH")

上面是用到插值的，比較慢，下面是取最近的整數的，比較快，但在某些極端的角度時效果很不好

import numpy as npimport netCDF4 as ncimport mathprint("IMPORT FINISH")##--read data--####--u--##da=nc.Dataset("4draw_average_u_3d")u=da.variables["u"][:]lon=da.variables["lon"][:]lat=da.variables["lat"][:]da.close()##--v--##da=nc.Dataset("4draw_average_v_3d")v=da.variables["v"][:]da.close()##--w--##da=nc.Dataset("4draw_average_w_3d")w=da.variables["w"][:]da.close()##--z--##da=nc.Dataset("4draw_average_z_3d")z=da.variables["z"][:]da.close()##--llj--##da=nc.Dataset("4draw_average_llj_3d")llj=da.variables["llj"][:]da.close()print("READ FINISH")##--line--##theta=0.2-np.pi/2#lat001=np.arange(50,160)#lon001=np.floor(lat001-100)*math.tan(theta)+165lon001=np.arange(80,170)lat001=np.floor((lon001-125)*0+130)#lat001=np.floor((lon001-125)/math.tan(theta)+130)lon001=list(lon001)lat001=list(lat001)if np.shape(lat001)[0]==np.shape(lon001)[0]:  ##--output--##  lon=np.array(lon)    lat=np.array(lat)  output_lon=lon[lon001,lat001]  output_lat=lat[lon001,lat001]  output_u=u[:,:,lon001,lat001]  output_v=v[:,:,lon001,lat001]  output_w=w[:,:,lon001,lat001]  output_z=z[:,:,lon001,lat001]  output_llj=llj[:,:,lon001,lat001]  ##--save data--##  da=nc.Dataset("4draw_crosssection_125_130","w",format="NETCDF4")  da.createDimension("heisize",np.shape(u)[1])  da.createDimension("xsize",np.shape(lat001)[0])  da.createDimension("timesize",np.shape(u)[0])  da.createDimension("onesize",1)  da.createVariable("u","f8",("timesize","heisize","xsize"))  da.createVariable("v","f8",("timesize","heisize","xsize"))  da.createVariable("w","f8",("timesize","heisize","xsize"))  da.createVariable("z","f8",("timesize","heisize","xsize"))  da.createVariable("llj","f8",("timesize","heisize","xsize"))  da.createVariable("lon","f8",("xsize",))  da.createVariable("lat","f8",("xsize",))  da.createVariable("lon001","f8",("xsize",))  da.createVariable("lat001","f8",("xsize",))  da.createVariable("theta","f8",("onesize",))  da.variables["u"][:]=output_u  da.variables["v"][:]=output_v  da.variables["w"][:]=output_w  da.variables["z"][:]=output_z  da.variables["llj"][:]=output_llj  da.variables["lat"][:]=output_lat  da.variables["lon"][:]=output_lon  da.variables["lat001"][:]=lat001  da.variables["lon001"][:]=lon001  da.variables["theta"][:]=theta  da.description="crosssection"  da.author="konghoiio"  da.createdate="2016-11-08"  da.close()else:  print("error happened")print("ALL FINISH")

以上的代碼可以得到crosssection的數據，為什么要把數據分析和畫圖的代碼分開呢？因為我是在服務器上做數據分析，在我的電腦上畫圖的ORZ

所以，下面是畫圖程序

def crosssection():  da=nc.Dataset("C:\\Users\\kong\\PycharmProjects\\southcnseapycode\\crosssection\\4draw_crosssection")  u=da.variables["u"][:]  v=da.variables["v"][:]  z=da.variables["z"][:]  w=da.variables["w"][:]  llj=da.variables["llj"][:]  lat=da.variables["lat"][:]  lon=da.variables["lon"][:]  lon001=da.variables["lon001"][:]  lat001=da.variables["lat001"][:]  #theta=da.variables["theta"][:]  theta=-np.pi/2  da.close()  print(z)  uv1=u*np.cos(theta)+v*np.sin(theta)##x  uv2=-u*np.sin(theta)+v*np.cos(theta)##y  lonlon=np.zeros((np.shape(z)[1],np.shape(lon)[0]))  latlat=np.zeros((np.shape(z)[1],np.shape(lat)[0]))  for i in range(0,np.shape(z)[1]):    lonlon[i,:]=lon    latlat[i,:]=lat  for i_hour in range(0,25):    print(i_hour)    plt.figure()    ##--contourf--##    cslev = np.arange(0,13)#13    cs=plt.contourf(latlat,z[i_hour,:,:],uv2[i_hour,:,:],cslev,cmap=cm.coolwarm,extend='both')    cs.cmap.set_over('k')    cs.cmap.set_under('b')    cbar = plt.colorbar(cs)    cbar.set_label('v\'(m/s)')    plt.ylim([0,4000])    plt.xlim([np.max(lat),np.min(lat)])    plt.xlabel("lat",size=20)    plt.ylabel("height(m)",size=20)    ##--vector--##    vectorsize=10    rho=3    Q = plt.quiver(latlat[:,::rho],z[i_hour,:,::rho],-uv1[i_hour,:,::rho],100000*w[i_hour,:,::rho],width=vectorsize,headwidth=vectorsize*2/3,headlength=vectorsize/2,headaxislength=vectorsize/2,minlength=0,units='dots',angles='xy', scale_units='xy', scale=50,minshaft=0.5)#100    qk = plt.quiverkey(Q, 0.1,1.1, 20, '20 m/s,uv&100000*w', labelpos='S')    ##--title--##    if i_hour==24:      plt.title("average",)    elif i_hour+7<24:      plt.title("%d UTC, %d localtime"%(i_hour,i_hour+7),fontsize="xx-large")    else:      plt.title("%d UTC, %d localtime"%(i_hour,i_hour+7-24),fontsize="xx-large")    plt.savefig('C:\\Users\\kong\\PycharmProjects\\southcnseapycode\\crosssection\\crosssection%d hour'%i_hour,dpi=500)    plt.close()

需要注意的是，橫座標是LON還是LAT呢？主要取決於你的截面與經度的夾角，夾角小於45度用LAT比較好，大於45用LON比較好。不過我的程序裏用LON比較方便，所以如非夾角很極端，我都建議用LON。不過上面的情況太極端了，所以用的是LAT

LON和LAT有什么區別嗎？

首先要注意

cs=plt.contourf(latlat,z[i_hour,:,:],uv2[i_hour,:,:],cslev,cmap=cm.coolwarm,extend='both')

plt.xlabel("lat",size=20)

 Q = plt.quiver(latlat[:,::rho],z[i_hour,:,::rho],-uv1[i_hour,:,::rho],100000*w[i_hour,:,::rho],width=vectorsize,headwidth=vectorsize*2/3,headlength=vectorsize/2,headaxislength=vectorsize/2,minlength=0,units='dots',angles='xy', scale_units='xy', scale=50,minshaft=0.5)#100

這裏面的latlat和"lat"要變。

另外，因於箭頭的正方向是指向低緯、高經度的。所以用LAT作橫坐標時，上面quiver中uv1前要有負號，為了使圖片右邊是正、左邊是負，要把坐標軸左右翻轉

plt.xlim([np.max(lat),np.min(lat)])

大概就是這樣了。

下面分別展示LAT、LON作橫坐標畫的圖。（所用的數據是不同的）