コンピュータを研究に使うための私的メモ。Python、Fortran、Ubuntu、etc...

極投影図、コンターのLocatorなど
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.ticker import *
import matplotlib.cm as cm
import matplotlib as mpl

def cplot(var,lon,lat):
    m = Basemap(projection='npstere',boundinglat=20,lon_0=180,\
            resolution='c',round=True)
    m.drawcoastlines(linewidth=0.5)
    m.drawmapboundary()
    m.drawmeridians(np.arange(0,360,30),linewidth=0.5,fontsize=14)
    m.drawparallels(np.arange(0,90,20),linewidth=0.5,fontsize=14)
    X,Y=m(lon,lat)
    m.contour(X,Y,var,colors='k',locator=MultipleLocator(100))

def mplot(var,ano,t,lon,lat):
    m = Basemap(projection='npstere',boundinglat=20,lon_0=180,\
            resolution='c',round=True)
    m.drawcoastlines(linewidth=0.5)
    m.drawmapboundary()
    m.drawmeridians(np.arange(0,360,30),linewidth=0.5,fontsize=14)
    m.drawparallels(np.arange(0,90,20),linewidth=0.5,fontsize=14)
    X,Y=m(lon,lat)
    if t*6<72:
        CM=m.contourf(X,Y,ano,np.arange(-80,81,10),cmap=cm.RdBu)
    elif t*6<24*5:
        CM=m.contourf(X,Y,ano,np.arange(-160,161,20),cmap=cm.RdBu)
    elif t*6<24*7:
        CM=m.contourf(X,Y,ano,np.arange(-200,201,25),cmap=cm.RdBu)
    else:
        CM=m.contourf(X,Y,ano,np.arange(-320,321,40),cmap=cm.RdBu)
    m.contour(X,Y,var,colors='k',locator=MultipleLocator(100))
    return CM

t=29
mpl.rc('figure.subplot',left=0.08,hspace=0.05,wspace=0,bottom=0.05,top=0.92)
plt.figure(figsize=(10,6))
plt.subplot(2,3,1)
mplot(hgt62[t,:,:],hgt62[t,:,:]-hgtgdas[t,:,:],t,lon,lat)
plt.text(0.8,0.9,'T62',fontsize=18,transform=plt.gca().transAxes)

plt.subplot(2,3,2)
mplot(hgt126[t,:,:],hgt126[t,:,:]-hgtgdas[t,:,:],t,lon,lat)
plt.text(0.8,0.9,'T126',fontsize=18,transform=plt.gca().transAxes)

plt.subplot(2,3,3)
mplot(hgt254[t,:,:],hgt254[t,:,:]-hgtgdas[t,:,:],t,lon,lat)
plt.text(0.8,0.9,'T254',fontsize=18,transform=plt.gca().transAxes)

plt.subplot(2,3,4)
CM=mplot(hgt574[t,:,:],hgt574[t,:,:]-hgtgdas[t,:,:],t,lon,lat)
plt.text(0.8,0.9,'T574',fontsize=18,transform=plt.gca().transAxes)
plt.text(0.6,-0.05,'contour interval = 100m',fontsize=12,transform=plt.gca().transAxes)
plt.colorbar(CM,plt.gcf().add_axes([0.9, 0.2, 0.015, 0.6]))
    
plt.suptitle('Z500 Forecast-Analysis  fhour='+str(t*6)+'h',fontsize=18)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from grads.ganum import GaNum

def rmse(err,lat):
    return np.sqrt(np.mean(np.average((err)**2,weights=np.cos(lat*np.pi/180.0),axis=1),axis=1))

fhour = np.arange(33)*6
plt.plot(fhour,rmse(hgt62-hgtgdas,lat), '-o',label='T62')
plt.plot(fhour,rmse(hgt126-hgtgdas,lat), '-o',label='T126')
plt.plot(fhour,rmse(hgt190-hgtgdas,lat), '-o',label='T190')
plt.plot(fhour,rmse(hgt254-hgtgdas,lat), '-o',label='T254')
plt.plot(fhour,rmse(hgt574-hgtgdas,lat), '-o',label='T574 (operational)')

from matplotlib.ticker import *
plt.gca().xaxis.set_minor_locator(MultipleLocator(6)) 
plt.gca().xaxis.set_major_locator(MultipleLocator(24)) 
plt.xlim([0,192])

plt.grid()
plt.legend(loc=2,numpoints=1)
plt.ylabel('Root Mean Square Error [m]')
plt.xlabel('Fourcast Hour [h]')
plt.title('20N-90N Z500 RMSE  2012/09/01 00UTC') 
plt.show()

情報

管理人/副管理人のみ編集できます