计算年分位数后绘画平均出错

lynpatty

大家好，有一个问题一直查不出错来，想要请教各位大神。
目前手头上有自己封装的NC文件，包含41年的tmax,tmin,precipitation的逐日降水资料及其经纬度。
现在是在每年每个季节的变量进行95的分位线计算（降水加一个对大于0的情况才进行95th分位），代码没有报错，但是画图出来令人匪夷所思。比如东亚的tmax，95th分位线上的夏季还是有-15度，然后降水的图出来也不是很好看，我想可能是哪个算法出了问题，还请各位大神帮我看看。
以下是代码：

1. import matplotlib.pyplot as plt
   
2. import matplotlib.cm as cm
   
3. import cartopy.crs as ccrs
   
4. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
   
5. import cartopy.feature as cfeature
   
6. import xarray as xr
   
7. import numpy as np
   
8. import pandas as pd
   
9. from pylab import *
   
10. from tqdm import tqdm
    
11. 
    
12. f = xr.open_dataset("/Volumes/Forests/era rusults/EA/data_day.nc")
    
13. ##f = xr.open_dataset('/Volumes/Forests/all sources data/1961-2018/CN05.1_Tmin_1961_2018_daily_025x025.nc',mask_and_scale=True)##我们输出的日数据###
    
14. tmin = f["tmin"] - 273.15
    
15. tmax = f["tmax"] - 273.15
    
16. pre = f["pre"] * 1000
    
17. 
    
18. lat = tmin.latitude
    
19. lon = tmin.longitude
    
20. def year_season_cal(var, year_s=1979,year_e=2019):
    
21.         year = np.arange(year_s,year_e+1,1)
    
22.         season_var = var[:4,:,:].expand_dims({"year":year},axis=0).copy()
    
23.         for i in range(0,len(year)):
    
24.                 ind, = np.where((var.time.dt.season== 'MAM')&(var.time.dt.year==year[i]))
    
25.                 season_var.values[i,0]= np.nanpercentile(var[ind],95,axis=0)
    
26.                 ind, = np.where((var.time.dt.season== 'JJA')&(var.time.dt.year==year[i]))
    
27.                 season_var.values[i,1]= np.nanpercentile(var[ind],95,axis=0)
    
28.                 ind, = np.where((var.time.dt.season== 'SON')&(var.time.dt.year==year[i]))
    
29.                 season_var.values[i,2] = np.nanpercentile(var[ind],95,axis=0)
    
30.                 ind = ind+91
    
31.                 season_var.values[i,3] = np.percentile(var[ind],95,axis=0)
    
32.                 return(season_var)
    
33. tmin_season = year_season_cal(tmin, year_s=1979,year_e=2019)
    
34. tmax_season = year_season_cal(tmax, year_s=1979,year_e=2019)
    
35. 
    
36. def pre_year_season_cal(var,year_s=1979,year_e=2019):
    
37.         year = np.arange(year_s,year_e+1,1)
    
38.         season_var = var[:4,:,:].expand_dims({"year":year},axis=0).copy()
    
39.         var_ = np.where(var==0,np.nan,var)
    
40.         for i in range(0,len(year)):
    
41.                 ind, = np.where((var.time.dt.season== 'MAM')&(var.time.dt.year==year[i]))
    
42.                 season_var.values[i,0] = np.nanpercentile(var_[ind],95,axis=0)
    
43.                 ind, = np.where((var.time.dt.season== 'JJA')&(var.time.dt.year==year[i]))
    
44.                 season_var.values[i,1] = np.nanpercentile(var_[ind],95,axis=0)
    
45.                 ind, = np.where((var.time.dt.season== 'SON')&(var.time.dt.year==year[i]))
    
46.                 season_var.values[i,2] = np.nanpercentile(var_[ind],95,axis=0)
    
47.                 ind = ind+91
    
48.                 season_var.values[i,3] = np.nanpercentile(var_[ind],95,axis=0)
    
49.                 return(season_var)
    
50. pre_season = pre_year_season_cal(pre, year_s=1979,year_e=2019)
    
51. 
    
52. tmin1 = np.zeros((4,len(lat),len(lon)))
    
53. for dd in range(4):
    
54.         tmin1[dd,:,:] = np.mean(tmax_season[:,dd,:,:],axis=0)
    
55. 
    
56. def plot(tmin1,lat,lon,cmap = plt.get_cmap("Spectral_r"),levels=19):
    
57.     fig, axes = plt.subplots(2, 2,figsize=(10,10),subplot_kw=dict(projection=ccrs.PlateCarree()))
    
58.     label = ["MAM", "JJA","SON","DJF"]
    
59.     for i, ax in enumerate(axes.flat):
    
60.         ax.set_extent([100, 145, 20, 50], crs=ccrs.PlateCarree())
    
61.         ax.set_xticks(np.linspace(100,145,5), crs=ccrs.PlateCarree())
    
62.         ax.set_yticks(np.linspace(20,50,5), crs=ccrs.PlateCarree())
    
63.         ax.add_feature(cfeature.LAND.with_scale('10m'))
    
64.         ax.add_feature(cfeature.OCEAN, zorder = 100, color = 'w')
    
65.         ax.add_feature(cfeature.COASTLINE.with_scale('10m'),lw=3)
    
66.         ax.add_feature(cfeature.BORDERS.with_scale('50m'),linestyle='-',lw=3)
    
67.         lon_formatter = LongitudeFormatter(number_format='.1f', degree_symbol=r'$^o
    
68. )
    
69.         lat_formatter = LatitudeFormatter(number_format='.1f', degree_symbol=r'$^o
    
70. )
    
71.         ax.xaxis.set_major_formatter(lon_formatter)
    
72.         ax.yaxis.set_major_formatter(lat_formatter)
    
73.         ax.tick_params(labelsize=8)
    
74.         cf = ax.contourf(lon, lat, tmin1[i], cmap=cmap,levels= levels, extend='both')
    
75.         ax.set_title(label[i], loc="left",fontsize=10)
    
76.         ax.set_title("Annual 95thTmax", loc="right",fontsize=10)
    
77.         cbar = plt.colorbar(cf, ax=ax, orientation='horizontal',
    
78.                         extendrect=False, pad=0.08, shrink=0.8)
    
79.         cbar.ax.tick_params(labelsize=8)
    
80.         plt.subplots_adjust(wspace=0.15, hspace=0.03)
    
81.     plt.savefig("Annual 95thTmax.png")
    
82. plot(tmin1,lat,lon)

复制代码

Lancelot

你可以用cdo处理一下，比较两者之间的差异，不过要注意numpy和cdo的分位数算法。祝好

lynpatty

Lancelot 发表于 2020-12-10 18:43
你可以用cdo处理一下，比较两者之间的差异，不过要注意numpy和cdo的分位数算法。祝好

我装了CDO但是没有用过，我想问你觉得这里我可以用CDO处理什么呢？

Lancelot

lynpatty 发表于 2020-12-10 18:49
我装了CDO但是没有用过，我想问你觉得这里我可以用CDO处理什么呢？

用cdo算分位数啊

lynpatty

Lancelot 发表于 2020-12-10 19:14
用cdo算分位数啊

可是还要对每年/每个季节划分呢？

lynpatty

Lancelot 发表于 2020-12-10 19:14
用cdo算分位数啊

hi你好，我昨天去查了一下CDO手册，找到了按年/季节划分然后计算分位数的办法，但是有一个问题能不能问一下你，手册说对于降水的划分也是要判断它是否大于0.1才进行计算，我想问这个对0.1的判断是CDO自动会进行的吗？还是我们要筛选一下导入进行计算呢？

Lancelot

lynpatty 发表于 2020-12-11 09:53
hi你好，我昨天去查了一下CDO手册，找到了按年/季节划分然后计算分位数的办法，但是有一个问题能不能问一 ...

应该不会，你可以用cdo 简单mask一下