PE paper 画图

CASEDIR='/Volumes/TU100Pro/1.个人研究工作专题/PEpaper2021DATA/'
from netCDF4 import Dataset
import numpy as np

import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
pltmap = plt.cm.rainbow
scale = '110m'
land = cfeature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',
facecolor=cfeature.COLORS['land'])

图2：画EN4融合数据的网格

Fig.2 TS profile locations (upper) and depths (below) of the period from 2010/01/01 to 2010/01/10.

NC = Dataset(CASEDIR+'EN_profile_thin_149399.nc')
lon = NC.variables['lon'][:]
lat = NC.variables['lat'][:]
depth = NC.variables['depth'][:]
temp = NC.variables['temp'][:]
salt= NC.variables['salt'][:]

rg = np.sum(temp<100,axis=1)
rgs = np.zeros(len(rg))
for j in range(len(rgs)):
    rgs[j] = depth[rg[j]-1]
# In[plt]
box = [0, 420, -80, 80]
scale = '50m'
xstep, ystep = 60, 40
fig = plt.figure(figsize=(10,12))
ax = fig.add_subplot(2, 1, 1, projection=ccrs.PlateCarree(central_longitude=180))
ax.set_extent(box, crs=ccrs.PlateCarree())
land = cfeature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',facecolor=cfeature.COLORS['land'])
ax.add_feature(land, facecolor='0.75')
ax.coastlines(scale)
ax.set_xticks(np.arange(box[0], box[1]+xstep,xstep), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(box[2], box[3]+ystep,ystep), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.scatter(lon,lat,s=3,color='k',transform=ccrs.PlateCarree(central_longitude=360))
ax.coastlines()
ax.set_global()
plt.title('Profile Locations 2010/1/1 - 2010/1/11',fontsize=16)
plt.xticks(fontsize=14);plt.yticks(fontsize=14)
ax = fig.add_subplot(4,1,3)
ax.bar(range(len(rgs)),rgs,color='k')
plt.yticks(np.arange(0,2500,500),fontsize=14)
plt.ylim(0,2200)
plt.xticks(np.arange(0,750,250),fontsize=14)
plt.xlim(0,500)
plt.grid(axis='y')
ax = plt.gca()
ax.invert_yaxis()
#plt.savefig(CASEDIR+'figure2.eps',format="eps")

图3: 敏感性实验结果

Fig.3 SST spread (left) and SSS spread (right) of the 1st, 4th, 7th,10th month from the parameter sensitivity experiment.

NC = Dataset(CASEDIR+'B_OSSE_sensitive.pop2009_std.nc')
LAT = NC.variables['TLAT'][:]
LON = NC.variables['TLONG'][:]
sst_std = NC.variables['SST'][:]
sss_std = NC.variables['SALT'][:,0,:,:]
MK = np.load(CASEDIR+'POP_region_mask.npz')
mk = MK['arr_0']

def interp_to_obs_grid(LON,LAT,SST,lon_obs,lat_obs):
    import numpy as np
#    SST is masked array
    [LONg,LATg] = np.meshgrid(lon_obs,lat_obs)
    points = np.array((LON.flatten(),LAT.flatten())).T
    value = SST.flatten();
    from scipy.interpolate import griddata
    SST_interped = griddata(points,value,(LONg,LATg),method = 'nearest')
    return SST_interped

lon = np.linspace(0.5,359.5,360)
lat = np.linspace(-79.5,80.5,160)
LONg,LATg = np.meshgrid(lon,lat)

SST_STD = np.zeros((12,160,360))
SSS_STD = np.zeros((12,160,360))

for j in range(12):
    sst_std[j][mk==0]=np.nan
    SST_STD[j] = interp_to_obs_grid(LON,LAT,sst_std[j],lon,lat)
    sss_std[j][mk==0]=np.nan
    SSS_STD[j] = interp_to_obs_grid(LON,LAT,sss_std[j],lon,lat)
    
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
pltmap = plt.cm.coolwarm
scale = '110m'
land = cfeature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',
                                    facecolor=cfeature.COLORS['land'])
SST_levels = np.arange(0,1.4,0.2)
SSS_levels = np.arange(0,0.6,0.1)
ylabels = ['1st Month','4th Month','7th Month','10th Month']
fig=plt.figure(figsize=(11,9))
for j in range(4):
    ax1 = plt.subplot(4,2,1+2*j,projection=ccrs.PlateCarree(central_longitude=180))
    ax1.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax1.add_feature(land, facecolor='0.75')
    ax1.coastlines(scale)
    if j==3:
        ax1.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    ax1.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax1.xaxis.set_major_formatter(lon_formatter)
    ax1.yaxis.set_major_formatter(lat_formatter)
    ax1.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SST_STD[j*3], levels = SST_levels,transform=ccrs.PlateCarree(),cmap = pltmap)
    plt.ylabel(ylabels[j],fontsize=14)
    if j==3:
        cbar = plt.colorbar(cor,cax = fig.add_axes([0.12, 0.05, 0.35, 0.02]),orientation='horizontal',pad=0.05,label='degC')
#    c=plt.contour(LONg,LATg,SST_STD[1+j*3],levels = SST_levels,transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
#    plt.clabel(c,inline=True,fontsize=10,fmt="%1.2f")
    
    if j==0:
        plt.title('SST spread',fontsize=16)

    ax2 = plt.subplot(4,2,2+2*j,projection=ccrs.PlateCarree(central_longitude=180))
    ax2.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax2.add_feature(land, facecolor='0.75')
    ax2.coastlines(scale)
    if j==3:
        ax2.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax2.xaxis.set_major_formatter(lon_formatter)
    ax2.yaxis.set_major_formatter(lat_formatter)
    ax2.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSS_STD[j*3],levels = SSS_levels, transform=ccrs.PlateCarree(),cmap = pltmap)
    if j==3:
        cbar = plt.colorbar(cor,cax = fig.add_axes([0.55, 0.05, 0.35, 0.02]),orientation='horizontal',pad=0.05,label='psu')
    if j==0:
        plt.title('SSS spread',fontsize=16)
#plt.savefig(CASEDIR+'figure3.eps',format='eps')

图4.状态估计稳定状态判定 Fig.4 RMSEs of prior and posterior innovations of the temperature and salinity for depths 10m,50m,200m, and 1000m.

fname= CASEDIR+'obs_diag_output.nc'
my_copy= 'rmse'
obstype1 = 'TEMPERATURE'
obstype2 = 'SALINITY'

from netCDF4 import Dataset
#import numpy as np

def byte_to_utf(byte_arr_in):
    output = []
    for j in range(len(byte_arr_in)):
        arrj = byte_arr_in[j]
        str_j = ''
        for a in arrj:
            str_j = str_j+a.decode('utf-8')
        output.append(str_j.rstrip())
    return output

DD = Dataset(fname)
COPIES = DD.variables['CopyMetaData'][:]
hlevel = DD.variables['hlevel'][:]
OBSTYPE = DD.variables['ObservationTypes'][:]
COPIES = byte_to_utf(COPIES)
OBSTYPE = byte_to_utf(OBSTYPE)
time = DD.variables['time'][:]

guess1 = DD.variables[obstype1+'_guess'][:]  # time,copy,level,region
analy1 = DD.variables[obstype1+'_analy'][:]
 
guess2 = DD.variables[obstype2+'_guess'][:] 
analy2 = DD.variables[obstype2+'_analy'][:]

TMPpri_sprd = guess1[1::,COPIES.index(my_copy),[0,1,2,3],0]
TMPpost_sprd = analy1[1::,COPIES.index(my_copy),[0,1,2,3],0]
SLTpri_sprd = guess2[1::,COPIES.index(my_copy),[0,1,2,3],0]
SLTpost_sprd = analy2[1::,COPIES.index(my_copy),[0,1,2,3],0]

time = time[1::]
TMPpri_sprd = TMPpri_sprd.reshape([-1,4])
TMPpost_sprd = TMPpost_sprd.reshape([-1,4])
SLTpri_sprd = 1000*SLTpri_sprd.reshape([-1,4])
SLTpost_sprd = 1000*SLTpost_sprd.reshape([-1,4])

# In[plt]
import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(12,9))
for p in range(4):
    plt.subplot(4,2,2*p+1)
    plt.plot(range(18), TMPpri_sprd[:,p],'b-o',lw=2,ms=8,label='prior')
    plt.plot(range(18), TMPpost_sprd[:,p],'r-s',lw=2,ms=8,label='post')
    plt.ylabel("%i"%hlevel[p]+'m',fontsize=12)
    if p==0:
        plt.title('Temperature(degC)',fontsize=14)
        plt.yticks(np.arange(0.3,1.7,0.3),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(5,20,5),[],fontsize=12);plt.xlim(-.3,17.3)
    elif p==1:
        plt.yticks(np.arange(0.5,1.5,0.2),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(5,20,5),[],fontsize=12);plt.xlim(-.3,17.3)
    elif p==2:
        plt.yticks(np.arange(0.45,0.7,0.05),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(5,20,5),[],fontsize=12);plt.xlim(-.3,17.3)
    else:
        plt.yticks(np.arange(0.2,0.4,0.05),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(5,20,5),['60','120','180'],fontsize=12)
        plt.xlabel('Days',fontsize=14);plt.xlim(-.3,17.3)
    plt.subplot(4,2,2*p+2)
    plt.plot(range(18), SLTpri_sprd[:,p],'b-o',lw=2,ms=8,label='prior')
    plt.plot(range(18), SLTpost_sprd[:,p],'r-s',lw=2,ms=8,label='post')
    if p==0:
        plt.title('Salinity(psu)',fontsize=14)
        plt.yticks(np.arange(0.10,0.26,0.03),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(5,20,5),[],fontsize=12);plt.xlim(-.3,17.3)
    elif p==1:
        plt.yticks(np.arange(0.08,0.15,0.02),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(5,20,5),[],fontsize=12);plt.xlim(-.3,17.3)
    elif p==2:
        plt.yticks(np.arange(0.055,0.08,0.005),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(0,20,5),[],fontsize=12);plt.xlim(-.3,17.3)
    else:
        plt.yticks(np.arange(0.023,0.03,0.002),fontsize=12)
        plt.grid();plt.legend(fontsize=12,ncol=2)
        plt.xticks(range(5,20,5),['60','120','180'],fontsize=12)
        plt.xlabel('Days',fontsize=14);plt.xlim(-.3,17.3)

#plt.savefig(CASEDIR+'figure4.eps',format='eps')

图5： RMSE Fig.5 RMSE of SST(upper),SSS(middle), and SSH(bottom) from control run (left) and state estimation (right) of the first 6 months.

import numpy as np
import os
from netCDF4 import Dataset
# ----- specify the case by names
this_year = 2009
ctl_case = "B_OSSE_control"
truth_case = "B_OSSE_truthobs"
se_case = "B_OSSE_SE"
pe_case = "B_OSSE_PE"
# ----- specify the dir of the nc file
truth_dir = "/Volumes/Seagate Expansion Drive/New_PE_expr/"+truth_case+'/'
ctl_dir = "/Volumes/Seagate Expansion Drive/New_PE_expr/"+ctl_case+'/'
se_dir="/Volumes/Seagate Expansion Drive/New_PE_expr/"+se_case+'/'
pe_dir="/Volumes/Seagate Expansion Drive/New_PE_expr/"+pe_case+'/'
# ------ read the netcdf data
TRUTH=Dataset(truth_dir+truth_case+'.pop'+str(this_year)+".nc") 
CTL = Dataset(ctl_dir+ctl_case+'.pop'+str(this_year)+".nc")
SE = Dataset(se_dir+se_case+'.pop'+str(this_year)+"_mean.nc")
PE = Dataset(pe_dir+pe_case+'.pop'+str(this_year)+"_mean.nc")
# ------ model grid
LON=TRUTH.variables["TLONG"][:]
LAT=TRUTH.variables["TLAT"][:]
depth = TRUTH.variables["z_t"][range(0,60,2)]/100 # every 2 levels
# ------- grid for evaluate and display
nmonth = 12; nlat = 160; nlon = 360;
lon = np.linspace(0.5,359.5,nlon)
lat = np.linspace(-79.5,79.5,nlat)
LONg,LATg = np.meshgrid(lon,lat)   # for plot
nlev=30
# function for interpolation
def interp_to_obs_grid(LON,LAT,SST,lon_obs,lat_obs):
    import numpy as np
#    SST is masked array
    [LONg,LATg] = np.meshgrid(lon_obs,lat_obs)
    points = np.array((LON.flatten(),LAT.flatten())).T
    value = SST.flatten();
    from scipy.interpolate import griddata
    SST_interped = griddata(points,value,(LONg,LATg),method = 'nearest')
    return SST_interped
# ----------function for compute errors
def computeErr(DATAdict,TRUTHdict):
    import numpy as np
    nmonth = 12; nlat = 160; nlon = 360; nlev = 30
    # interpolate to uniform grids
    lon = np.linspace(0.5,359.5,nlon)
    lat = np.linspace(-79.5,79.5,nlat)
    LON=TRUTHdict.variables["TLONG"][:]
    LAT=TRUTHdict.variables["TLAT"][:]
    TEMP0 = TRUTHdict.variables["TEMP"][:,range(0,60,2)]
    TEMP1 = DATAdict.variables["TEMP"][:,range(0,60,2)]
    SALT0 = TRUTHdict.variables["SALT"][:,range(0,60,2)]
    SALT1 = DATAdict.variables["SALT"][:,range(0,60,2)]
    SST0 = TRUTHdict.variables["SST"][:]
    SST1 = DATAdict.variables["SST"][:]
    SSH0 = TRUTHdict.variables["SSH"][:]
    SSH1 = DATAdict.variables["SSH"][:]
    TEMPerr = np.zeros((nmonth,nlev,nlat,nlon))
    SALTerr = np.zeros((nmonth,nlev,nlat,nlon))
    SSTerr = np.zeros((nmonth,nlat,nlon))
    SSHerr = np.zeros((nmonth,nlat,nlon))
    for n in range(nmonth):
        for l in range(nlev):
            TEMPerr[n,l] = interp_to_obs_grid(LON,LAT,TEMP1[n,l]-TEMP0[n,l],lon,lat)
            SALTerr[n,l] = interp_to_obs_grid(LON,LAT,SALT1[n,l]-SALT0[n,l],lon,lat)
        SSTerr[n] = interp_to_obs_grid(LON,LAT,SST1[n]-SST0[n],lon,lat)
        SSHerr[n] = interp_to_obs_grid(LON,LAT,SSH1[n]-SSH0[n],lon,lat)
        print('procedeed month'+str(n+1))
    return TEMPerr,SALTerr,SSTerr,SSHerr
# In[load data]: compute or load
# ------for control
if not os.path.exists(ctl_dir + ctl_case+str(this_year)+"data.npz"):
    TEMPerrC,SALTerrC,SSTerrC,SSHerrC = computeErr(CTL,TRUTH)
    np.savez(ctl_dir + ctl_case+str(this_year)+"data.npz",TEMPerrC,SALTerrC,SSTerrC,SSHerrC)
else:
    LDATA = np.load(ctl_dir + ctl_case+str(this_year)+"data.npz")
    TEMPerrC = LDATA['arr_0'];SALTerrC = LDATA['arr_1'];SSTerrC = LDATA['arr_2'];SSHerrC = LDATA['arr_3']
TEMPerrC[TEMPerrC>100]=np.nan; SALTerrC[SALTerrC>100]=np.nan; 
SSTerrC[SSTerrC>100]=np.nan;SSHerrC[SSHerrC>100]=np.nan; 
# ------for SE
if not os.path.exists(se_dir + se_case+str(this_year)+"data.npz"):
    TEMPerrS,SALTerrS,SSTerrS,SSHerrS = computeErr(SE,TRUTH)
    np.savez(se_dir + se_case+str(this_year)+"data.npz",TEMPerrS,SALTerrS,SSTerrS,SSHerrS)
else:
    LDATA = np.load(se_dir + se_case+str(this_year)+"data.npz")
    TEMPerrS = LDATA['arr_0'];SALTerrS = LDATA['arr_1'];SSTerrS = LDATA['arr_2'];SSHerrS = LDATA['arr_3']
TEMPerrS[TEMPerrS>100]=np.nan; SALTerrS[SALTerrS>100]=np.nan; 
SSTerrS[SSTerrS>100]=np.nan;SSHerrS[SSHerrS>100]=np.nan; 
# ------for PE
if not os.path.exists(pe_dir + pe_case+str(this_year)+"data.npz"):
    TEMPerrP,SALTerrP,SSTerrP,SSHerrP = computeErr(PE,TRUTH)
    np.savez(pe_dir + pe_case+str(this_year)+"data.npz",TEMPerrP,SALTerrP,SSTerrP,SSHerrP)
else:
    LDATA = np.load(pe_dir + pe_case+str(this_year)+"data.npz")
    TEMPerrP = LDATA['arr_0'];SALTerrP = LDATA['arr_1'];SSTerrP = LDATA['arr_2'];SSHerrP = LDATA['arr_3']
TEMPerrP[TEMPerrP>100]=np.nan; SALTerrP[SALTerrP>100]=np.nan; 
SSTerrP[SSTerrP>100]=np.nan;SSHerrP[SSHerrP>100]=np.nan;
# In[plt surface]
import matplotlib.pyplot as plt
import numpy as np

if True:
    SSTlevels = np.arange(0,3,0.5)
    SSSlevels = np.arange(0,1,0.2)
    SSHlevels = np.arange(0,18,3)
    
    SSTrmseC = np.sqrt(np.mean(np.square(SSTerrC[range(2,6)]),axis=0))
    SSTrmseS = np.sqrt(np.mean(np.square(SSTerrS[range(2,6)]),axis=0))
    SSTrmseP = np.sqrt(np.mean(np.square(SSTerrP[range(2,6)]),axis=0))
    
    SSHrmseC = np.sqrt(np.mean(np.square(SSHerrC[range(2,6)]),axis=0))
    SSHrmseS = np.sqrt(np.mean(np.square(SSHerrS[range(2,6)]),axis=0))
    SSHrmseP = np.sqrt(np.mean(np.square(SSHerrP[range(2,6)]),axis=0))
    
    SSSerrC = SALTerrC[:,0];
    SSSerrS = SALTerrS[:,0];SSSerrP = SALTerrP[:,0]
    SSSrmseC = np.sqrt(np.mean(np.square(SSSerrC[range(6)]),axis=0))
    SSSrmseS = np.sqrt(np.mean(np.square(SSSerrS[range(6)]),axis=0))
    SSSrmseP = np.sqrt(np.mean(np.square(SSSerrP[range(6)]),axis=0))
    ##
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
    pltmap = plt.cm.coolwarm
    scale = '110m'
    land = cfeature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',
                                    facecolor=cfeature.COLORS['land'])

    fig=plt.figure(figsize=(12,8))
    ax1 = plt.subplot(3,2,1,projection=ccrs.PlateCarree(central_longitude=180))
    ax1.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax1.add_feature(land, facecolor='0.75')
    ax1.coastlines(scale)
#    ax1.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    ax1.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax1.xaxis.set_major_formatter(lon_formatter)
    ax1.yaxis.set_major_formatter(lat_formatter)
    ax1.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSTrmseC, levels=SSTlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
#    cbar = plt.colorbar(cor,orientation='vertical',shrink=0.7,pad=0.05,label='degC')
#    c=plt.contour(LONg,LATg,SSTrmseC,levels = SSTlevels,transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
#    plt.clabel(c,inline=True,fontsize=10,fmt="%1.2f")
    plt.title('Control Run',fontsize=16)
    plt.ylabel('SST',fontsize=15)
    
    ax2 = plt.subplot(3,2,2,projection=ccrs.PlateCarree(central_longitude=180))
    ax2.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax2.add_feature(land, facecolor='0.75')
    ax2.coastlines(scale)
#    ax2.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
#    ax2.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
#    lat_formatter = LatitudeFormatter()
    ax2.xaxis.set_major_formatter(lon_formatter)
    ax2.yaxis.set_major_formatter(lat_formatter)
    ax2.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSTrmseS, levels=SSTlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
    plt.title('Assimilated',fontsize=16)
    cbar = plt.colorbar(cor,cax = fig.add_axes([0.93, 0.67, 0.01, 0.2]),orientation='vertical',label='degC')
#    c=plt.contour(LONg,LATg,SSTrmseS,levels = SSTlevels,transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
#    plt.clabel(c,inline=True,fontsize=10,fmt="%1.2f")
    
    
    ax3 = plt.subplot(3,2,3,projection=ccrs.PlateCarree(central_longitude=180))
    ax3.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax3.add_feature(land, facecolor='0.75')
    ax3.coastlines(scale)
#    ax3.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    ax3.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax3.xaxis.set_major_formatter(lon_formatter)
    ax3.yaxis.set_major_formatter(lat_formatter)
    ax3.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSSrmseC, levels=SSSlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
#    cbar = plt.colorbar(cor,orientation='vertical',shrink=0.7,pad=0.05,label='psu')
#    c=plt.contour(LONg,LATg,SSSrmseC,levels = SSSlevels,transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
#    plt.clabel(c,inline=True,fontsize=10,fmt="%1.2f")
    plt.ylabel('SSS',fontsize=15)
    
    ax4 = plt.subplot(3,2,4,projection=ccrs.PlateCarree(central_longitude=180))
    ax4.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax4.add_feature(land, facecolor='0.75')
    ax4.coastlines(scale)
#    ax4.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
#    ax4.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax4.xaxis.set_major_formatter(lon_formatter)
    ax4.yaxis.set_major_formatter(lat_formatter)
    ax4.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSSrmseS, levels=SSSlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
#    cbar = plt.colorbar(cor,orientation='vertical',shrink=0.7,pad=0.05,label='psu')
#    c=plt.contour(LONg,LATg,SSSrmseS,levels = SSSlevels,transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
#    plt.clabel(c,inline=True,fontsize=10,fmt="%1.2f")
    cbar = plt.colorbar(cor,cax = fig.add_axes([0.93, 0.4, 0.01, 0.2]),orientation='vertical',label='psu')
    
    ax5 = plt.subplot(3,2,5,projection=ccrs.PlateCarree(central_longitude=180))
    ax5.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax5.add_feature(land, facecolor='0.75')
    ax5.coastlines(scale)
    ax5.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    ax5.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax5.xaxis.set_major_formatter(lon_formatter)
    ax5.yaxis.set_major_formatter(lat_formatter)
    ax5.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSHrmseC,levels=SSHlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
#    cbar = plt.colorbar(cor,orientation='vertical',shrink=0.7,pad=0.05,label='cm')
#    c=plt.contour(LONg,LATg,SSHrmseC,levels = SSHlevels,transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
#    plt.clabel(c,inline=True,fontsize=10,fmt="%1.2f")
    plt.ylabel('SSH',fontsize=15)
    
    ax6 = plt.subplot(3,2,6,projection=ccrs.PlateCarree(central_longitude=180))
    ax6.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax6.add_feature(land, facecolor='0.75')
    ax6.coastlines(scale)
    ax6.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
#    ax6.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax6.xaxis.set_major_formatter(lon_formatter)
    ax6.yaxis.set_major_formatter(lat_formatter)
    ax6.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSHrmseS, levels=SSHlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
#    cbar = plt.colorbar(cor,orientation='vertical',shrink=0.7,pad=0.05,label='cm')
#    c=plt.contour(LONg,LATg,SSHrmseS,levels = SSHlevels,transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
#    plt.clabel(c,inline=True,fontsize=10,fmt="%1.2f")
    cbar = plt.colorbar(cor,cax = fig.add_axes([0.93, 0.13, 0.01, 0.2]),orientation='vertical',label='cm')

#plt.savefig(CASEDIR+'figure5.eps',format='eps')

图6： Fig.6 (a).the ratio of posterior error standard deviation to prior error standard deviation of the parameter ensemble; (b-d).the posterior ensemble member 1 - 3; for the first parameter estimation step.

def interp_to_obs_grid(LON,LAT,SST,lon_obs,lat_obs):
    import numpy as np
#    SST is masked array
    [LONg,LATg] = np.meshgrid(lon_obs,lat_obs)
    points = np.array((LON.flatten(),LAT.flatten())).T
    value = SST.flatten();
    from scipy.interpolate import griddata
    SST_interped = griddata(points,value,(LONg,LATg),method = 'nearest')
    return SST_interped

lon = np.linspace(0.5,359.5,360)
lat = np.linspace(-79.5,80.5,160)
LONg,LATg = np.meshgrid(lon,lat)

#%%
NC = Dataset(CASEDIR+"postassim_sd.nc")
VPARA = NC.variables['VMIX_PARA'][:]
LON = NC.variables['TLON'][:]
LAT = NC.variables['TLAT'][:]
NC2 = Dataset(CASEDIR+"preassim_sd.nc")
VPARA0 = NC2.variables['VMIX_PARA'][:]
ratio = interp_to_obs_grid(LON, LAT, VPARA/VPARA0, lon, lat)
paras = np.zeros([3,160,360])
for j in range(3):
    NCj = Dataset(CASEDIR+'postassim_member_000'+str(j+1)+'.nc')
    vtmp = NCj.variables['VMIX_PARA'][:]
    paras[j] = interp_to_obs_grid(LON, LAT, vtmp, lon, lat)

pltmap=plt.cm.coolwarm_r
#%%
fig=plt.figure(figsize=(12,5))
ax1 = plt.subplot(2,2,1,projection=ccrs.PlateCarree(central_longitude=180))
ax1.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
ax1.add_feature(land, facecolor='0.75')
ax1.coastlines(scale)
#ax1.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
ax1.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
ax1.coastlines(resolution=scale)
plt.title(r'$\sigma^u/\sigma^p$',fontsize=12)
cor = plt.contourf(LONg,LATg,ratio,levels=np.arange(0.48,1.04,0.05), transform=ccrs.PlateCarree(),cmap = pltmap)
cbar = plt.colorbar(shrink=0.5)
c=plt.contour(LONg,LATg,ratio,levels = np.array([0.83]),transform=ccrs.PlateCarree(), colors='black',linestyles='solid')
plt.clabel(c,inline=True,fontsize=11,fmt="%1.2f")

plt.tight_layout(h_pad=0)
lvl_para=np.arange(0.16,0.24,0.015)
for j in range(3):
    ax2 = plt.subplot(2,2,j+2,projection=ccrs.PlateCarree(central_longitude=180))
    ax2.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax2.add_feature(land, facecolor='0.75')
    ax2.coastlines(scale)
    if j>0:
        ax2.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    if j==1:
        ax2.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax2.xaxis.set_major_formatter(lon_formatter)
    ax2.yaxis.set_major_formatter(lat_formatter)
    ax2.coastlines(resolution=scale)
    plt.title(r'$\theta_'+str(j+1)+'$',fontsize=12)
    cor = plt.contourf(LONg,LATg,paras[j],levels=lvl_para, transform=ccrs.PlateCarree(),cmap = plt.cm.rainbow)
    cbar = plt.colorbar(shrink=0.5)
    c=plt.contour(LONg,LATg,ratio,levels = np.array([0.83]),transform=ccrs.PlateCarree(), colors='black',linestyles='solid')
    plt.clabel(c,inline=True,fontsize=11,fmt="%1.2f")
#plt.savefig(CASEDIR+'figure6.eps',format='eps')

图7 Fig.7 The spatial varying inflation values for temperature (left) and salinity (right) in 5m (a,b), 100m (c,d) and 200m (e,f) depth; (g) the inflation values for the two-dimensional parameter field for the first analysis step.

def interp_to_obs_grid(LON,LAT,SST,lon_obs,lat_obs):
    import numpy as np
#    SST is masked array
    [LONg,LATg] = np.meshgrid(lon_obs,lat_obs)
    points = np.array((LON.flatten(),LAT.flatten())).T
    value = SST.flatten();
    from scipy.interpolate import griddata
    SST_interped = griddata(points,value,(LONg,LATg),method = 'nearest')
    return SST_interped

lon = np.linspace(0.5,359.5,360)
lat = np.linspace(-79.5,80.5,160)
LONg,LATg = np.meshgrid(lon,lat)

#%%
NC = Dataset(CASEDIR+"postassim_priorinf_mean.nc")

SST_inf = NC.variables['TEMP_CUR'][[0,9,19]]
SSS_inf = NC.variables['SALT_CUR'][[0,9,19]]
PARA_inf = NC.variables['VMIX_PARA'][:] 
LON = NC.variables['TLON'][:]
LAT = NC.variables['TLAT'][:]

SST_INF = np.zeros((3,160,360))
SSS_INF = np.zeros((3,160,360))
for j in range(3):
    SST_INF[j] = interp_to_obs_grid(LON,LAT,SST_inf[j],lon,lat)
    SSS_INF[j] = interp_to_obs_grid(LON,LAT,SSS_inf[j],lon,lat)
PARA_INF = interp_to_obs_grid(LON, LAT,PARA_inf , lon, lat)
#%%
SST_levels = np.arange(1,3.3,0.4)
SSS_levels = np.arange(1,3.3,0.4)
ylabels = ['5m','100m','200m']
fig=plt.figure(figsize=(11,9))
for j in range(3):
    ax1 = plt.subplot(4,2,1+2*j,projection=ccrs.PlateCarree(central_longitude=180))
    ax1.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax1.add_feature(land, facecolor='0.75')
    ax1.coastlines(scale)
    if j==2:
        ax1.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    ax1.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax1.xaxis.set_major_formatter(lon_formatter)
    ax1.yaxis.set_major_formatter(lat_formatter)
    ax1.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SST_INF[j], levels = SST_levels,transform=ccrs.PlateCarree(),cmap = plt.cm.Reds)
    if j==0:
        plt.title('TEMPERATURE',fontsize=14)
    plt.ylabel(ylabels[j],fontsize=14)
    if j==2:
        cbar = plt.colorbar(cor,cax = fig.add_axes([0.49, 0.37,  0.01, 0.45,]),orientation='vertical',pad=0.05)

    
    ax2 = plt.subplot(4,2,2+2*j,projection=ccrs.PlateCarree(central_longitude=180))
    ax2.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
    ax2.add_feature(land, facecolor='0.75')
    ax2.coastlines(scale)
    if j==2:
        ax2.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
    lon_formatter = LongitudeFormatter(zero_direction_label=False)
    lat_formatter = LatitudeFormatter()
    ax2.xaxis.set_major_formatter(lon_formatter)
    ax2.yaxis.set_major_formatter(lat_formatter)
    ax2.coastlines(resolution=scale)
    cor = plt.contourf(LONg,LATg,SSS_INF[j],levels = SSS_levels, transform=ccrs.PlateCarree(),cmap = plt.cm.Blues)
    if j==0:
        plt.title('SALINITY',fontsize=14)
    if j==2:
        cbar = plt.colorbar(cor,cax = fig.add_axes([0.91, 0.37,  0.01, 0.45,]),orientation='vertical',pad=0.05)
    
ax3 = plt.subplot(4,1,4,projection=ccrs.PlateCarree(central_longitude=180))
ax3.set_extent([0,420,-70,70], crs=ccrs.PlateCarree())
ax3.add_feature(land, facecolor='0.75')
ax3.coastlines(scale)
ax3.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
ax3.set_yticks([-60, -30, 0, 30, 60], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax3.xaxis.set_major_formatter(lon_formatter)
ax3.yaxis.set_major_formatter(lat_formatter)
ax3.coastlines(resolution=scale)
cor = plt.contourf(LONg,LATg,PARA_INF,levels = SSS_levels, transform=ccrs.PlateCarree(),cmap = plt.cm.Greens)
cbar = plt.colorbar(cor,shrink=0.8,label='inflation')
plt.ylabel('Para. field',fontsize=14)
#plt.savefig(CASEDIR+'figure7.eps',format='eps')

Text(0, 0.5, 'Para. field')

图8：

f = open(CASEDIR+"PEtxt/parameters1.0.txt")
lines = f.readlines()
ListP = lines[0::4]
N = len(ListP)
for j in range(N):
    ListP[j] = ListP[j]+lines[4*j+1]+lines[4*j+2]+lines[4*j+3]
    ListP[j] = ListP[j][1:-2]
import numpy as np
Data1 = np.zeros((N,20))
for n in range(N):
    tmp = ListP[n].split()
    for k in range(20):
        Data1[n,k] = float(tmp[k])
# data assimilation starts at the 10th cycle
Data = np.zeros((N+17,len(Data1[0])))
for j in range(N+17):
    if j<17:
        Data[j] = Data1[0];
    else:
        Data[j] = Data1[j-17];
Dmean = np.mean(Data,axis=1)
Dsd = np.std(Data,axis=1)

import matplotlib.pyplot as plt
plt.figure(figsize=(12,6))
plt.subplot(211)
plt.plot(Dmean,'r-o',lw=2,ms=8,label='Para. Ens. Mean')
plt.plot(0.16*np.ones(N+17),'k',lw=3,label='Truth')
plt.plot(0.10*np.ones(N+17),'k--',lw=3,label='Initial Guess')
plt.plot(Dmean-Dsd,'r:',lw=2)
plt.plot(Dmean+Dsd,'r:',lw=2)
plt.xticks(range(0,N+17,9),[],fontsize=15);plt.xlim(-0.5,36.5)
plt.yticks(np.arange(0.04,0.23,0.06),fontsize=15)
plt.grid()
plt.legend(fontsize=14,ncol=3,loc=9)
plt.title('Parameter Ensemble Evolution',fontsize=15)

plt.subplot(212)
plt.plot(0.16-Dmean,'r-o',lw=3,label='Absolute Error')
plt.plot(Dsd,'g-s',lw=3,label='Ensemble Spread')
plt.xticks(range(8,N+17,9),['Apr','Jul','Oct','Jan'],fontsize=15);plt.xlim(-0.5,36.5)
plt.yticks(np.arange(0,0.1,0.02),fontsize=15)
plt.grid()
plt.legend(fontsize=14,ncol=3,loc=8)
plt.title('Ensemble Spread and Error',fontsize=15)

#plt.savefig(CASEDIR+'figure8.eps',format='eps')

Text(0.5, 1.0, 'Ensemble Spread and Error')

图9：

import numpy as np
from netCDF4 import Dataset
import datetime as dtm

d0 = dtm.datetime(2009,6,30)
N = 18
lvl_idx = [0,5,20,39]

spdT0 = np.zeros([18,40]); spdS0 = np.zeros([18,40]);
spdT1 = np.zeros([18,40]); spdS1 = np.zeros([18,40]);

DT= np.load(CASEDIR+'/spds.npz')
spdT0 = DT['arr_0']
spdT1 = DT['arr_1']
spdS0 = DT['arr_2']*1000
spdS1 = DT['arr_3']*1000
rateT = np.roll(spdT1,-1,axis=0)/spdT0
rateS = np.roll(spdS1,-1,axis=0)/spdS0
time = np.arange(1,19)
zt = DT['arr_4']

import matplotlib.pyplot as plt
plt.figure(figsize=(12,9))
for p in range(3):
    ax = plt.subplot(4,2,2*p+1)
    ax2 = ax.twinx()
    for j in range(17):
        ax.plot([time[j],time[j+1]],[spdT0[j,lvl_idx[p]],spdT1[j+1,lvl_idx[p]]],'ro--',ms=5)
    ax2.bar(time[1:18]-0.5,rateT[0:17,lvl_idx[p]],fc='k',width=0.5)
    if p==0:
        ax.set_xlim(0.5,18.5)
        ax.set_title('Gobal mean temperature spread',fontsize=15)
        ax.set_ylabel('5m',fontsize=14)
        ax.set_ylim(0,0.28);ax.set_yticks([0.15,0.2,0.25]);ax.grid(axis='y',linestyle='-',fillstyle='left')
        ax.set_yticklabels([0.15,0.20,0.25],fontsize=14)
        ax2.set_ylim(1,4);ax2.set_yticks([1,1.5,2]);ax2.grid(axis='y',linestyle='--',fillstyle='right')
        ax2.set_yticklabels([1.0,1.5,2.0],fontsize=14)
        ax.set_xlabel('Time',fontsize=14)
        plt.xticks([1,6,11,16],['06-30','08-19','10-08','11-27'],fontsize=15)
    elif p==1:
        ax.set_xlim(0.5,18.5)
        ax.set_ylabel('50m',fontsize=14)
        ax.set_ylim(0,0.23);ax.set_yticks([0.1,0.15,0.2]);ax.grid(axis='y',linestyle='-',fillstyle='left')
        ax.set_yticklabels([0.1,0.15,0.2],fontsize=14)
        ax2.set_ylim(1,4);ax2.set_yticks([1,1.5,2]);ax2.grid(axis='y',linestyle='--',fillstyle='right')
        ax2.set_yticklabels([1.0,1.5,2.0],fontsize=14)
        ax.set_xlabel('Time',fontsize=14)
        plt.xticks([1,6,11,16],['06-30','08-19','10-08','11-27'],fontsize=15)
#    elif p==2:
    else:
        ax.set_xlim(0.5,18.5)
        ax.set_ylabel('200m',fontsize=14)
        ax.set_ylim(0,0.1);ax.set_yticks([0.05,0.07,0.09]);ax.grid(axis='y',linestyle='-',fillstyle='left')
        ax.set_yticklabels([0.05,0.07,0.09],fontsize=14)
        ax2.set_ylim(1,4);ax2.set_yticks([1,1.5,2]);ax2.grid(axis='y',linestyle='--',fillstyle='right')
        ax2.set_yticklabels([1.0,1.5,2.0],fontsize=14)
        ax.set_xlabel('Time',fontsize=14)
        plt.xticks([1,6,11,16],['06-30','08-19','10-08','11-27'],fontsize=15)

    ax = plt.subplot(4,2,2*p+2)
    ax2 = ax.twinx()
    for j in range(17):   
        ax.plot([time[j],time[j+1]],[spdS0[j,lvl_idx[p]],spdS1[j+1,lvl_idx[p]]],'bo--',ms=5)
    ax2.bar(time[1:18]-0.5,rateS[0:17,lvl_idx[p]],fc='k',width=0.5) 
    if p==0:
        ax.set_xlim(0.5,18.5)
        ax.set_title('Gobal mean salinity spread',fontsize=15)
        ax.set_ylim(0,0.1);ax.set_yticks([0.04,0.06,0.08]);ax.grid(axis='y',linestyle='-',fillstyle='left')
        ax.set_yticklabels([0.04,0.06,0.08],fontsize=14)
        ax2.set_ylim(1,2);ax2.set_yticks([1,1.2,1.4]);ax2.grid(axis='y',linestyle='--',fillstyle='right')
        ax2.set_yticklabels([1,1.2,1.4],fontsize=14)
        ax2.set_ylabel('Ratio',fontsize=14)
        ax.set_xlabel('Time',fontsize=14)
        plt.xticks([1,6,11,16],['06-30','08-19','10-08','11-27'],fontsize=15)
    elif p==1:
        ax.set_xlim(0.5,18.5)
        ax.set_ylim(0,0.04);ax.set_yticks([0.02,0.03,0.04]);ax.grid(axis='y',linestyle='-',fillstyle='left')
        ax.set_yticklabels([0.02,0.03,0.04],fontsize=14)
        ax2.set_ylim(1,2);ax2.set_yticks([1,1.2,1.4]);ax2.grid(axis='y',linestyle='--',fillstyle='right')
        ax2.set_yticklabels([1,1.2,1.4],fontsize=14)
        ax2.set_ylabel('Ratio',fontsize=14)
        ax.set_xlabel('Time',fontsize=14)
        plt.xticks([1,6,11,16],['06-30','08-19','10-08','11-27'],fontsize=15)
#    elif p==2:
    else:
        ax.set_xlim(0.5,18.5)
        ax.set_ylim(0,0.013);ax.set_yticks([0.007,0.01,0.013]);ax.grid(axis='y',linestyle='-',fillstyle='left')
        ax.set_yticklabels(['.007','.010','.013'],fontsize=14)
        ax2.set_ylim(1,2);ax2.set_yticks([1,1.2,1.4]);ax2.grid(axis='y',linestyle='--',fillstyle='right')
        ax2.set_yticklabels([1,1.2,1.4],fontsize=14)
        ax2.set_ylabel('Ratio',fontsize=14)
        ax.set_xlabel('Time',fontsize=14)
        plt.xticks([1,6,11,16],['06-30','08-19','10-08','11-27'],fontsize=15)
    plt.tight_layout()
    plt.subplot(4,2,7)
    plt.plot(zt[range(40)],rateT[0],'k')
    plt.xlim(0,1000);plt.xticks(np.arange(0,1250,250),fontsize=14)
    plt.xlabel('Depth',fontsize=14)
    plt.ylabel('Ratio',fontsize=14)
    plt.yticks([1.0,1.2,1.4,1.6,1.8],fontsize=14)
    plt.grid(axis='y',linestyle='--',fillstyle='left')
    plt.title('Temperature spread rate of 07-10',fontsize=15)
    plt.subplot(4,2,8)
    plt.plot(zt[range(40)],rateS[0],'k')
    plt.xlim(0,1000);plt.xticks(np.arange(0,1250,250),fontsize=14)
    plt.yticks([1.0,1.1,1.2,1.3,1.4],fontsize=14)
    plt.xlabel('Depth',fontsize=14)
    plt.title('Salinity spread rate of 07-10',fontsize=15)
    plt.grid(axis='y',linestyle='--',fillstyle='left')
    
#plt.savefig(CASEDIR+'figure9.eps',format='eps')

<ipython-input-9-733e9d499b89>:95: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.subplot(4,2,7)
<ipython-input-9-733e9d499b89>:103: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.subplot(4,2,8)
<ipython-input-9-733e9d499b89>:95: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.subplot(4,2,7)
<ipython-input-9-733e9d499b89>:103: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  plt.subplot(4,2,8)

图10：

PARAS_mean = []
PARAs_sd = []
# read the data and convert to arrays
f = open(CASEDIR+"parameters_inf1.1bnd0.016.txt")
lines = f.readlines()
ListP = lines[0::4]
N = len(ListP)
for j in range(N):
    ListP[j] = ListP[j]+lines[4*j+1]+lines[4*j+2]+lines[4*j+3]
    ListP[j] = ListP[j][1:-2]
import numpy as np
Data1 = np.zeros((N,20))
for n in range(N):
    tmp = ListP[n].split()
    for k in range(20):
        Data1[n,k] = float(tmp[k])
# data assimilation starts at the 10th cycle
Data = np.zeros((N+17,len(Data1[0])))
for j in range(N+17):
    if j<17:
        Data[j] = Data1[0];
    else:
        Data[j] = Data1[j-17];
Dmean = np.mean(Data,axis=1)
Dsd = np.std(Data,axis=1)
PARAS_mean.append(Dmean)
PARAs_sd.append(Dsd)
f = open(CASEDIR+"parameters_inf1.2bnd0.016.txt")
lines = f.readlines()
ListP = lines[0::4]
N = len(ListP)
for j in range(N):
    ListP[j] = ListP[j]+lines[4*j+1]+lines[4*j+2]+lines[4*j+3]
    ListP[j] = ListP[j][1:-2]
import numpy as np
Data1 = np.zeros((N,20))
for n in range(N):
    tmp = ListP[n].split()
    for k in range(20):
        Data1[n,k] = float(tmp[k])
# data assimilation starts at the 10th cycle
Data = np.zeros((N+17,len(Data1[0])))
for j in range(N+17):
    if j<17:
        Data[j] = Data1[0];
    else:
        Data[j] = Data1[j-17];
Dmean = np.mean(Data,axis=1)
Dsd = np.std(Data,axis=1)
PARAS_mean.append(Dmean)
PARAs_sd.append(Dsd)
f = open(CASEDIR+"parameters_inf1.3bnd0.016.txt")
lines = f.readlines()
ListP = lines[0::4]
N = len(ListP)
for j in range(N):
    ListP[j] = ListP[j]+lines[4*j+1]+lines[4*j+2]+lines[4*j+3]
    ListP[j] = ListP[j][1:-2]
import numpy as np
Data1 = np.zeros((N,20))
for n in range(N):
    tmp = ListP[n].split()
    for k in range(20):
        Data1[n,k] = float(tmp[k])
# data assimilation starts at the 10th cycle
Data = np.zeros((N+17,len(Data1[0])))
for j in range(N+17):
    if j<17:
        Data[j] = Data1[0];
    else:
        Data[j] = Data1[j-17];
Dmean = np.mean(Data,axis=1)
Dsd = np.std(Data,axis=1)
PARAS_mean.append(Dmean)
PARAs_sd.append(Dsd)
f = open(CASEDIR+"parameters.txt")
lines = f.readlines()
ListP = lines[0::4]
N = len(ListP)
for j in range(N):
    ListP[j] = ListP[j]+lines[4*j+1]+lines[4*j+2]+lines[4*j+3]
    ListP[j] = ListP[j][1:-2]
import numpy as np
Data1 = np.zeros((N,20))
for n in range(N):
    tmp = ListP[n].split()
    for k in range(20):
        Data1[n,k] = float(tmp[k])
# data assimilation starts at the 10th cycle
Data = np.zeros((N+17,len(Data1[0])))
for j in range(N+17):
    if j<17:
        Data[j] = Data1[0];
    else:
        Data[j] = Data1[j-17];
Dmean = np.mean(Data,axis=1)
Dsd = np.std(Data,axis=1)
PARAS_mean.append(Dmean)
PARAs_sd.append(Dsd)
# In[plt]: plot figures  
inf_alpha = [1.1,1.2,1.3,1.5]     
import matplotlib.pyplot as plt
plt.figure(figsize=(10,10))
for j in range(4):
    plt.subplot(4,1,j+1)
    plt.plot(PARAS_mean[j],'r',lw=3,label='Para. Ens. Mean')
    plt.plot(PARAS_mean[j]+PARAs_sd[j],'r:',lw=3)
    plt.plot(PARAS_mean[j]-PARAs_sd[j],'r:',lw=3)
    
    plt.plot(0.16*np.ones(N+17),'k',lw=4,label='Truth')
    plt.plot(0.10*np.ones(N+17),'k--',lw=4,label='State Est. Guess')
    plt.xticks(range(8,N+17,9),[],fontsize=14)
    plt.yticks([0.04,0.1,0.16,0.22],fontsize=14)
    
    plt.xlim(0,N+17)
    plt.grid()
    plt.title(r'$\alpha=$ '+str(inf_alpha[j]),fontsize=15)
    if j == 0:
        plt.legend(fontsize=14,ncol=3)
    if j == 3:
        plt.xticks(np.arange(8,N+17,9),['Apr','Jul','Oct','Jan','Apr','Jul','Oct','Jan'],fontsize=14)

#plt.savefig(CASEDIR+'figure10.eps',format='eps')

import numpy as np
import os
from netCDF4 import Dataset
# ----- specify the case by names
this_year = 2010
truth_case = "B_OSSE_truthobs"
se_case = "B_OSSE_SE"
pe_case = "B_OSSE_PE"
# ----- specify the dir of the nc file
truth_dir = "/Volumes/Seagate Expansion Drive/New_PE_expr/"+truth_case+'/'
se_dir="/Volumes/Seagate Expansion Drive/New_PE_expr/"+se_case+'/'
pe_dir="/Volumes/Seagate Expansion Drive/New_PE_expr/"+pe_case+'_inf1.2bnd0.016/'
# ------ read the netcdf data
TRUTH=Dataset(truth_dir+truth_case+'.pop'+str(this_year)+".nc") 
SE = Dataset(se_dir+se_case+'.pop'+str(this_year)+"_mean.nc")
PE = Dataset(pe_dir+pe_case+'.pop'+str(this_year)+"_mean.nc")
# ------ model grid
LON=TRUTH.variables["TLONG"][:]
LAT=TRUTH.variables["TLAT"][:]
depth = TRUTH.variables["z_t"][range(0,60,2)]/100 # every 2 levels
# ------- grid for evaluate and display
nmonth = 12; nlat = 160; nlon = 360;
lon = np.linspace(0.5,359.5,nlon)
lat = np.linspace(-79.5,79.5,nlat)
LONg,LATg = np.meshgrid(lon,lat)   # for plot
nlev=30
# function for interpolation
def interp_to_obs_grid(LON,LAT,SST,lon_obs,lat_obs):
    import numpy as np
#    SST is masked array
    [LONg,LATg] = np.meshgrid(lon_obs,lat_obs)
    points = np.array((LON.flatten(),LAT.flatten())).T
    value = SST.flatten();
    from scipy.interpolate import griddata
    SST_interped = griddata(points,value,(LONg,LATg),method = 'nearest')
    return SST_interped
# ----------function for compute errors
def computeErr(DATAdict,TRUTHdict):
    import numpy as np
    nmonth = 12; nlat = 160; nlon = 360; nlev = 30
    # interpolate to uniform grids
    lon = np.linspace(0.5,359.5,nlon)
    lat = np.linspace(-79.5,79.5,nlat)
    LON=TRUTHdict.variables["TLONG"][:]
    LAT=TRUTHdict.variables["TLAT"][:]
    TEMP0 = TRUTHdict.variables["TEMP"][:,range(0,60,2)]
    TEMP1 = DATAdict.variables["TEMP"][:,range(0,60,2)]
    SALT0 = TRUTHdict.variables["SALT"][:,range(0,60,2)]
    SALT1 = DATAdict.variables["SALT"][:,range(0,60,2)]
    SST0 = TRUTHdict.variables["SST"][:]
    SST1 = DATAdict.variables["SST"][:]
    SSH0 = TRUTHdict.variables["SSH"][:]
    SSH1 = DATAdict.variables["SSH"][:]
    TEMPerr = np.zeros((nmonth,nlev,nlat,nlon))
    SALTerr = np.zeros((nmonth,nlev,nlat,nlon))
    SSTerr = np.zeros((nmonth,nlat,nlon))
    SSHerr = np.zeros((nmonth,nlat,nlon))
    for n in range(nmonth):
        for l in range(nlev):
            TEMPerr[n,l] = interp_to_obs_grid(LON,LAT,TEMP1[n,l]-TEMP0[n,l],lon,lat)
            SALTerr[n,l] = interp_to_obs_grid(LON,LAT,SALT1[n,l]-SALT0[n,l],lon,lat)
        SSTerr[n] = interp_to_obs_grid(LON,LAT,SST1[n]-SST0[n],lon,lat)
        SSHerr[n] = interp_to_obs_grid(LON,LAT,SSH1[n]-SSH0[n],lon,lat)
        print('procedeed month'+str(n+1))
    return TEMPerr,SALTerr,SSTerr,SSHerr

if not os.path.exists(se_dir + se_case+str(this_year)+"data.npz"):
    TEMPerrS,SALTerrS,SSTerrS,SSHerrS = computeErr(SE,TRUTH)
    np.savez(se_dir + se_case+str(this_year)+"data.npz",TEMPerrS,SALTerrS,SSTerrS,SSHerrS)
else:
    LDATA = np.load(se_dir + se_case+str(this_year)+"data.npz")
    TEMPerrS = LDATA['arr_0'];SALTerrS = LDATA['arr_1'];SSTerrS = LDATA['arr_2'];SSHerrS = LDATA['arr_3']
TEMPerrS[TEMPerrS>100]=np.nan; SALTerrS[SALTerrS>100]=np.nan; 
SSTerrS[SSTerrS>100]=np.nan;SSHerrS[SSHerrS>100]=np.nan; 
# ------for PE
if not os.path.exists(pe_dir + pe_case+str(this_year)+"data.npz"):
    TEMPerrP,SALTerrP,SSTerrP,SSHerrP = computeErr(PE,TRUTH)
    np.savez(pe_dir + pe_case+str(this_year)+"data.npz",TEMPerrP,SALTerrP,SSTerrP,SSHerrP)
else:
    LDATA = np.load(pe_dir + pe_case+str(this_year)+"data.npz")
    TEMPerrP = LDATA['arr_0'];SALTerrP = LDATA['arr_1'];SSTerrP = LDATA['arr_2'];SSHerrP = LDATA['arr_3']
TEMPerrP[TEMPerrP>100]=np.nan; SALTerrP[SALTerrP>100]=np.nan; 
SSTerrP[SSTerrP>100]=np.nan;SSHerrP[SSHerrP>100]=np.nan;
# In[plt surface]
tRMSES = np.zeros((nmonth,nlev,2));tRMSEP = np.zeros((nmonth,nlev,2));

Tlatrange=range(80,140);Slatrange=range(50,110)

for n in range(nmonth):
    for l in range(nlev):
        tRMSES[n,l,0] = np.sqrt(np.nanmean(np.square(TEMPerrS[n,l,Tlatrange])))
        tRMSEP[n,l,0] = np.sqrt(np.nanmean(np.square(TEMPerrP[n,l,Tlatrange])))
        
        tRMSES[n,l,1] = np.sqrt(np.nanmean(np.square(SALTerrS[n,l,Slatrange])))
        tRMSEP[n,l,1] = np.sqrt(np.nanmean(np.square(SALTerrP[n,l,Slatrange])))

plt.figure(figsize=(15,9))
plt.tight_layout()
for j in range(6):
    plt.subplot(2,6,j+1)
    plt.plot(tRMSES[j+6,:,0],depth,'k--',lw=3,label='State Est.')
    plt.plot(tRMSEP[j+6,:,0],depth,'r',lw=3,label='Para. Est.')
    plt.ylim(0,400);plt.xlim(0,0.5)
    ax = plt.gca()
    ax.invert_yaxis()
    plt.title(str(j+7)+'th month',fontsize=14)
    plt.grid()
    if j==0:
        plt.yticks(np.arange(0,400,100),fontsize=12)
        plt.legend(fontsize=12,loc=8)
        plt.ylabel('TEMP RMSE', fontsize=14)
    else:
        plt.yticks(np.arange(0,400,100),[],fontsize=12)    
    plt.xticks([0,0.2,0.4],fontsize=12)
    if j==2:
        plt.xlabel('degC', fontsize=12)
        
    plt.subplot(2,6,j+7)
    plt.plot(tRMSES[j+6,:,1],depth,'k--',lw=3,label='State Est.')
    plt.plot(tRMSEP[j+6,:,1],depth,'r',lw=3,label='Para. Est.')
    plt.ylim(0,400);plt.xlim(0,0.13)
    if j==0:
        plt.yticks(np.arange(0,400,100),fontsize=12)
        plt.legend(fontsize=12,loc=8)
        plt.ylabel('SALT RMSE', fontsize=14)
    else:
        plt.yticks(np.arange(0,400,100),[],fontsize=12)   
    plt.xticks([0,0.05,0.1],fontsize=12)
    ax = plt.gca()
    ax.invert_yaxis()
    plt.grid()
    if j==2:
        plt.xlabel('psu', fontsize=12)
#plt.savefig(CASEDIR+'figure11.eps',format='eps')

# SSTlevels = np.arange(0,1.6,0.2)
# SSSlevels = np.arange(0,0.25,0.04)
# SSHlevels = np.arange(0,18,3)
# fig=plt.figure(figsize=(12,10))
# ax1 = plt.subplot(4,1,1,projection=ccrs.PlateCarree(central_longitude=180))
# ax1.set_extent([120,240,35,60], crs=ccrs.PlateCarree())
# ax1.add_feature(land, facecolor='0.75')
# ax1.coastlines(scale)
# ax1.set_xticks([ 120, 160, 200, 240], crs=ccrs.PlateCarree())
# ax1.set_yticks([40,50,60], crs=ccrs.PlateCarree())
# lon_formatter = LongitudeFormatter(zero_direction_label=False)
# lat_formatter = LatitudeFormatter()
# ax1.xaxis.set_major_formatter(lon_formatter)
# ax1.yaxis.set_major_formatter(lat_formatter)
# ax1.coastlines(resolution=scale)
# cor = plt.contourf(LONg,LATg,abs(SSTerrS[8]), levels=SSTlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
# c=plt.contour(LONg,LATg,abs(SSTerrS[8]), levels=np.arange(0,1.6,0.4),transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
# plt.ylabel('State Est.',fontsize=15)

# ax1 = plt.subplot(4,1,2,projection=ccrs.PlateCarree(central_longitude=180))
# ax1.set_extent([120,240,35,60], crs=ccrs.PlateCarree())
# ax1.add_feature(land, facecolor='0.75')
# ax1.coastlines(scale)
# ax1.set_xticks([ 120, 160, 200, 240], crs=ccrs.PlateCarree())
# ax1.set_yticks([40,50,60], crs=ccrs.PlateCarree())
# lon_formatter = LongitudeFormatter(zero_direction_label=False)
# lat_formatter = LatitudeFormatter()
# ax1.xaxis.set_major_formatter(lon_formatter)
# ax1.yaxis.set_major_formatter(lat_formatter)
# ax1.coastlines(resolution=scale)
# cor = plt.contourf(LONg,LATg,abs(SSTerrP[8]), levels=SSTlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
# c=plt.contour(LONg,LATg,abs(SSTerrP[8]), levels=np.arange(0,1.6,0.4),transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
# plt.ylabel('Para. Est.',fontsize=15)
# cbar = plt.colorbar(cor,cax = fig.add_axes([0.9, 0.55,  0.01, 0.3,]),orientation='vertical',pad=0.05,label='degC')

# ax2 = plt.subplot(4,1,3,projection=ccrs.PlateCarree(central_longitude=180))
# ax2.set_extent([140,260,-10,15], crs=ccrs.PlateCarree())
# ax2.add_feature(land, facecolor='0.75')
# ax2.coastlines(scale)
# ax2.set_xticks([140, 180, 220, 260], crs=ccrs.PlateCarree())
# ax2.set_yticks([-10,0,10], crs=ccrs.PlateCarree())
# lon_formatter = LongitudeFormatter(zero_direction_label=False)
# lat_formatter = LatitudeFormatter()
# ax2.xaxis.set_major_formatter(lon_formatter)
# ax2.yaxis.set_major_formatter(lat_formatter)
# ax2.coastlines(resolution=scale)
# cor = plt.contourf(LONg,LATg,abs(SSSerrS[6]), levels=SSSlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
# c=plt.contour(LONg,LATg,abs(SSSerrS[6]), levels=np.arange(0.04,0.25,0.08),transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
# plt.ylabel('State Est.',fontsize=15)

# ax2 = plt.subplot(4,1,4,projection=ccrs.PlateCarree(central_longitude=180))
# ax2.set_extent([140,260,-10,15], crs=ccrs.PlateCarree())
# ax2.add_feature(land, facecolor='0.75')
# ax2.coastlines(scale)
# ax2.set_xticks([140, 180, 220, 260], crs=ccrs.PlateCarree())
# ax2.set_yticks([-10,0,10], crs=ccrs.PlateCarree())
# lon_formatter = LongitudeFormatter(zero_direction_label=False)
# lat_formatter = LatitudeFormatter()
# ax2.xaxis.set_major_formatter(lon_formatter)
# ax2.yaxis.set_major_formatter(lat_formatter)
# ax2.coastlines(resolution=scale)
# cor = plt.contourf(LONg,LATg,abs(SSSerrP[6]), levels=SSSlevels, transform=ccrs.PlateCarree(),cmap = pltmap)
# c=plt.contour(LONg,LATg,abs(SSSerrP[6]), levels=np.arange(0.04,0.25,0.08),transform=ccrs.PlateCarree(), colors='black',linestyles='dashed')
# plt.ylabel('Para. Est.',fontsize=15)
# cbar = plt.colorbar(cor,cax = fig.add_axes([0.9, 0.15,  0.01, 0.3,]),orientation='vertical',pad=0.05,label='psu')
# #plt.savefig(CASEDIR+'figure12.eps',format='eps')