任务：
将连续5天日平均气温均大于20°的一段时期识别为热浪,请根据ERA5的温度数据(t2mdaily2000_2009.nc)
计算并画图:
()请根据热浪定义,识别热浪事件,画出区域内平均的每年热浪发生天数的时间序列图
(2)计算并画出该区域内多年平均的热浪发生频率的空间分布图
6月1日晚上七点
# 读取nc文件并对其中的数据进行处理
# 作者： 林子轩
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
file = "t2m_daily_2000_2009.nc"
dataset = nc.Dataset(file)
print(dataset.variables.keys())
dict_keys(['time', 'longitude', 'latitude', 't2m'])
dataset.variables['time']
<class 'netCDF4._netCDF4.Variable'>
int64 time(time)
units: days since 2000-01-01 00:00:00
calendar: proleptic_gregorian
unlimited dimensions:
current shape = (3653,)
filling on, default _FillValue of -9223372036854775806 used
思路：[-2,-1,0,1,2] 都大于
dataset.variables['t2m'][:].data.shape # (3653, 41, 41)
(3653, 41, 41)
F = dataset.variables['t2m'][:].data
C = (F-273.15)
C.shape # (3653, 41, 41)
(3653, 41, 41)
"""获取热浪事件数 边界处理 """
def get_hot_event(daily_mean):
hot_day = []
for day , daily_temp in enumerate(daily_mean):
flag = 1
if day-2 >0 and day+3 <len(daily_mean): # 边界处理
for value in daily_mean[day-2 :day + 3]:
if value>20:
pass
else: flag = 0
# 连续5天平均气温大于 20
if flag: hot_day.append(day)
hot_event = []
pre_day = -1
# 找出value 的间隔数目
for today in hot_day:
if pre_day+1 != today:
hot_event.append(today)
pre_day = today
return np.array(hot_event)
"""记录每个小格子的热浪事件"""
hot_event = []
m = 0
for i in range(C.shape[1]):
for j in range(C.shape[2]):
hot_event.append(get_hot_event(C[:,i,j]))
len(hot_event)
1681
C[0,0,0]
-11.762939
i = 0
for day in hot_event[0]:
if i<10:
print(day , C[day-3:day+4,0,0])
else:
i +=1
163 [19.088257 21.79602 22.580292 23.014648 25.512573 20.593933 20.535492]
188 [17.784729 20.894104 23.290924 20.544159 21.878235 23.6734 24.788239]
205 [18.21811 21.270294 21.804504 22.227905 26.531586 27.927368 20.94983 ]
214 [18.639557 22.12616 20.301636 21.015503 23.227081 21.94403 21.292206]
226 [19.437744 20.89215 21.815125 20.136353 21.203613 20.159393 20.445374]
505 [16.130768 23.219543 23.087463 25.246887 26.518097 24.569977 13.239929]
519 [18.883636 21.406586 24.294617 25.607208 26.684052 26.392822 22.641754]
539 [19.452667 20.308105 20.803162 23.898315 26.18634 27.986511 23.315186]
546 [19.401062 21.483887 22.973969 21.589722 20.883453 22.64447 23.634277]
567 [18.37326 21.88556 23.170593 22.402771 27.270538 21.742554 17.227844]
574 [16.449036 21.1279 22.202423 22.693329 21.60086 23.513794 23.81897 ]
587 [19.512177 21.209808 23.659363 22.463745 23.507538 23.100403 16.152496]
882 [19.63861 23.739594 24.222015 24.397888 27.190216 27.090515 26.345947]
896 [17.962952 20.662537 20.779755 25.017883 24.5896 26.819244 20.154358]
909 [19.34314 22.636627 23.046478 23.306488 24.160187 21.143005 22.62793 ]
939 [19.14972 21.904663 22.014618 23.202972 20.428986 24.262512 23.873718]
956 [18.26941 21.31607 22.247986 21.822083 23.94049 23.286377 23.893799]
1263 [18.926208 20.341125 23.651031 23.141907 24.649841 22.807007 24.097717]
1287 [18.308838 22.392517 23.902283 23.10962 24.355957 22.429504 20.05072 ]
1299 [19.494873 20.832947 21.797913 22.190735 20.83731 24.325745 18.059387]
1305 [18.059387 22.30899 23.939514 20.76883 21.660767 21.603333 21.469666]
1622 [18.225891 20.652771 22.849396 24.801208 24.111084 23.168854 15.556946]
1639 [19.217041 22.070831 24.870422 24.904846 23.761322 20.550842 22.43335 ]
1649 [19.123535 20.562714 22.628479 24.058258 25.182617 20.21759 14.65387 ]
1657 [18.335693 20.508545 21.148132 23.440857 23.740662 23.713135 22.488708]
1994 [19.066101 22.98529 24.35376 23.367676 26.14038 23.769836 25.995056]
2012 [19.335907 22.342987 22.361511 25.87857 27.335144 25.336792 22.585327]
2019 [19.692017 21.05957 23.328308 27.485107 25.21466 26.175903 27.940918]
2030 [17.798218 20.183594 22.77069 24.826477 26.013916 26.520752 16.441925]
2037 [19.191498 21.164795 21.96875 21.825806 20.3609 21.898163 18.655457]
2043 [18.655457 22.196777 22.234833 23.762451 23.5607 22.76297 23.20111 ]
2363 [19.80655 21.684113 23.613281 26.316833 21.887756 20.814667 19.925842]
2373 [17.513672 21.375732 22.631226 25.641418 25.685211 25.683014 25.051514]
2392 [18.964417 21.771362 22.594849 23.904327 22.627258 21.446686 19.26651 ]
2406 [19.786499 20.418335 21.8266 24.046814 24.432373 24.383057 24.112762]
2713 [11.905518 20.007751 24.65091 25.70398 25.573822 27.24408 27.678802]
2733 [16.869232 22.500793 24.602783 26.461273 23.897339 20.789978 21.544495]
2749 [19.807526 21.71228 23.192749 25.789124 21.985046 22.168427 25.280823]
2764 [19.343018 20.968567 23.313751 22.76709 20.79596 21.012512 19.687073]
2774 [19.230713 20.642029 22.738586 22.450134 22.087006 23.390778 21.835876]
2797 [19.909271 22.236145 20.437958 21.407166 21.310883 20.434448 18.272156]
3081 [15.416443 21.70221 24.323975 26.37146 24.309204 25.457642 22.990936]
3099 [16.111816 20.548737 22.0177 20.483673 20.66095 20.951843 19.131073]
3105 [19.131073 20.441711 23.585144 26.00711 24.963104 20.424683 21.167908]
3120 [16.94873 20.329987 22.922974 21.73581 22.335114 22.567535 21.709167]
3137 [17.005005 22.501709 24.763855 25.965271 26.148926 26.188477 25.24646 ]
3454 [19.465668 21.71231 25.943329 23.034332 26.98529 20.277588 18.487183]
3470 [19.935974 20.07077 21.951233 24.943634 26.963562 20.711304 20.484833]
3490 [17.653564 24.35376 27.473145 22.63086 20.67273 20.049988 23.285248]
3497 [19.476807 21.67038 22.468048 21.43155 22.830872 22.543304 21.201294]
3511 [18.357788 21.246521 24.802368 25.484863 26.013428 21.275818 19.89148 ]
"""获取每个小格子的热浪事件个数"""
hot_num = []
for value in hot_event: # 1681
hot_num.append(len(value))
hot_num = np.array(hot_num).reshape(41,41)
hot_num[0,0]
51
plt.imshow(hot_num)
<matplotlib.image.AxesImage at 0x282f668bca0>
def get_figure(zz):
x = dataset.variables['longitude'][:].data
y = dataset.variables['latitude'][:].data
xx,yy = np.meshgrid(x, y)
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=16) #导入本地汉字字库
plt.subplots(figsize=(12,8))#该方法会返回画图对象 和坐标对象ax
plt.contourf(xx, yy, zz,FontProperties =font)
plt.title("Hot Events Freq",fontsize=16)
plt.ylabel("Latitude(deg)",fontsize=16)
plt.xlabel("Longitude(deg)",fontsize=16)
plt.colorbar()
plt.savefig(f'figures/hot_event{zz.shape}')
plt.gca().yaxis.set_major_formatter(mticker.FormatStrFormatter('%.0f°N'))
plt.gca().xaxis.set_major_formatter(mticker.FormatStrFormatter('%.0f°E'))
plt.show()
格点平均 每 365天
"""获取每年的每个小格子个数"""
per_year_hot = np.zeros((10,1681))
for m,high in enumerate(range(365,len(C),365)):
low = high - 365
for j,value in enumerate(hot_event): # 每一个小格子 切成10份\ 每个小格子位置和对应数据
flag = (np.array(low<value)+0) * (np.array(value<high)+0)
temp = flag*value
num = len(np.unique(temp))-1
per_year_hot[m,j] = num
mean_figure2 = np.mean(per_year_hot,axis=1)
per_year_hot = per_year_hot.reshape(10,41,41)
mean_year = np.mean(per_year_hot,axis=0)
mean_year.shape,mean_figure2.shape
((41, 41), (10,))
mean_figure2,1
(array([3.96371208, 3.70136823, 3.92325996, 4.37775134, 4.62522308,
3.76264128, 4.06662701, 4.51338489, 3.99286139, 4.84533016]),
1)
get_figure(mean_year)
C:\Users\11958\AppData\Local\Temp\ipykernel_29960\498370392.py:8: UserWarning: The following kwargs were not used by contour: 'FontProperties'
plt.contourf(xx, yy, zz,FontProperties =font)
def get_figure2(freq):
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=16) #导入本地汉字字库
plt.subplots(figsize=(12,8))#该方法会返回画图对象 和坐标对象ax
plt.title("Frequency of Hot Events From 2000 To 2009",fontsize=16)
plt.ylabel("Frequency(times/year)",fontsize=16)
plt.xlabel("year",fontsize=16)
plt.plot(list(range(2000,2010)),freq,'r--')
plt.plot(list(range(2000,2010)),freq,'ro')
plt.savefig(f'figures/hot_event{freq.shape}')
plt.show()
get_figure2(mean_figure2)
1101
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
