【Python气象绘图临摹】处理数据(上):读入输出nc数据、截取夏季/冬季数据、ButterWorth带通滤波、计算方差
前言
2022.9學(xué)習(xí)繪圖
利用python進(jìn)行氣象繪圖,本文為學(xué)習(xí)繪制期間記錄筆記,分為上、下兩部分:處理數(shù)據(jù)和圖像繪制。處理數(shù)據(jù)流程:讀入olr資料,截取夏季/冬季數(shù)據(jù),進(jìn)行10-30dButterWorth帶通濾波,計(jì)算方差,輸出nc資料。
繪制用到資料orl.day.mean.nc
參考文獻(xiàn)
[1] Qian Y , Hsu P C , Kazuyoshi K . New real-time indices for the quasi-biweekly oscillation over the Asian summer monsoon region[J]. Climate Dynamics, 2019.
->圖1 Variance of 10–30-day BP-filtered OLR (shading; units:W2m-4) and propagation vectors of 10–30-day-filtered OLR for a the annual-mean state, b boreal summer (May– Oct.) and c boreal winter (Nov.– Apr.) during 1980–2012.
代碼
濾波原文用的是30年數(shù)據(jù),考慮到內(nèi)存和計(jì)算時(shí)長本文僅用了1999年1年的數(shù)據(jù)計(jì)算,大概分布也是一致的,如果感興趣也可以用30年算算,會更趨近于原文。
import time time_start = time.process_time() import numpy as np import xarray as xr import datetime as dt #處理datetime64格式數(shù)據(jù)使用,可以解析時(shí)間格式 from scipy import signal import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import cartopy.mpl.ticker as cticker ####################################################### ## 讀取olr數(shù)據(jù) ####################################################### f = xr.open_dataset('E:/LACS/learn/python/replot/data/olr.day.mean.nc') # print(f) #xarray.Dataset 對象 twStrt, twLast = 1999,1999 LatS, LatN, LonL, LonR= -45.0, 45.0, 0.0, 360.0 olr = f['olr'].loc['%s-01-01'%twStrt:'%s-12-31'%twLast].loc[:,LatS:LatN,LonL:LonR] # olr.shape>>>(365, 21, 360)####################################################### ## 10–30-day BP-filter ####################################################### # 去線性趨勢;設(shè)置為'constant',則為去平均值,即為求距平 olr_1030 = signal.detrend(olr, axis=0, type='linear')#設(shè)置為'constant',則為去平均值,即為求距平# ButterWorth帶通濾波 b, a = signal.butter(6, [1/30.,1/10.], 'bandpass') olr_1030 = signal.filtfilt(b, a, olr_1030,axis = 0) ####################################################### ## 計(jì)算方差Varizances ####################################################### olr_1030 = xr.DataArray(data=olr_1030, #轉(zhuǎn)為數(shù)據(jù)數(shù)組DataArraydims=["time", "lat", "lon"],coords=[olr.time,olr.lat,olr.lon],attrs=dict(description="olr_10-30filtered",units="W2/m4")) ## 全年 ## var_all = np.var(olr_1030,axis=0) ## 夏季5-10月 ## # 索引1999-1999年之間5-10月 經(jīng)度0-360 緯度0-20N的數(shù)據(jù) olr_sum = olr_1030.loc[olr_1030.time.dt.month.isin([5,6,7,8,9,10])] #.shape>>>(184, 21, 360) var_sum = np.var(olr_sum,axis=0) # ## 冬季11-4月 ## # # 索引1999-1999年之間11-次年4月 經(jīng)度0-360 緯度0-20N的數(shù)據(jù) olr_win = olr_1030.loc[olr_1030.time.dt.month.isin([11,12,1,2,3,4])]#.shape>>>(182 184, 181, 360) var_win = np.var(olr_win,axis=0) print(var_sum) print(olr)####################################################### ## 輸出成nc文件 ####################################################### ds = xr.Dataset({'var_all':var_all, 'var_sum':var_sum, 'var_win':var_win})# 將數(shù)據(jù)數(shù)組轉(zhuǎn)為數(shù)據(jù)集 ds.to_netcdf('./var_all_sum_win.nc', mode='w') # 使用相對路徑也可以time_end = time.process_time() print('\nRunning time: %s Seconds'%(time_end-time_start))總結(jié)
以上是生活随笔為你收集整理的【Python气象绘图临摹】处理数据(上):读入输出nc数据、截取夏季/冬季数据、ButterWorth带通滤波、计算方差的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 淘宝图片搜索API / item_sea
- 下一篇: python读取cad_python3读