我在 python 方面没有经验,我正在将 NCL 脚本转换为 python,希望 python 运行得更快。环顾四周,我没有找到我认为是 NCL 脚本中最简单计算的答案。看看如何完成更艰难的计算,我也没有找到如何在 python 中完成这些的答案。
大部分计算是在将 3 维变量转换为 1 维变量并查询它们在数组空间中的值和位置之后完成的。知道了 t 个变量在数组空间中的位置,我们可以得到对应于 t 个变量整数值的 p 个变量值。
计算如下:
- 将 p 变量中的值设置为其默认值 _FillValue,
- 计算网格点的数量(体积) t 变量中每个可能的整数值都会出现一个值(时间和空间的总和),
- 计算 t 变量中每个可能的整数值的开始和结束时间索引,
- 将持续时间计算为 t 变量中结束时间和开始时间之间的差值(+1,因为数字),
- 计算 t 变量中每个可能的整数值的平均(时空)纬度和经度,
- 计算 t 变量中每个可能的整数值的面积(体积/持续时间),
- 从 p 变量计算平均 p,其中它在空间时间上对应于 t 变量中每个可能的整数值,并且
- 从 p 变量计算 p 百分位数,其中它在时空中对应于 t 变量中的每个可能的整数值。
所有这些计算都将值保存在一维数组中,维度大小等于 t 变量中的最大整数值。例如,at 变量可能有 0 到 100 之间的整数。0 整数值被忽略,因此示例中的每个一维数组应该有 100 个值;(100 卷、100 个开始时间、100 个结束时间等)。
最后,将所有一维数组写入(制表符分隔的)文本文件,每一列都是一维数组。
;===================================================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;===================================================================
begin
;===============================================================
begTime = get_cpu_time()
; Data I/O and data names
; T output, and raw data input
f_t = addfile("t_in.nc","r")
f_p = addfile("p_in.nc","r")
; Data variables
time = f_t->time
lat = f_t->lat
lon = f_t->lon
t_var = f_t->t
p_var = f_p->p
p_fix = p_var
p_fix = where(p_var.eq.9.96921e+36, p_var@_FillValue, p_var)
delete(p_var)
p_var = p_fix
delete(p_fix)
; t = 0, is not measured
; Compute volume, start and end time indices, delta time, centroid lat and centroid lon, area, and percentiles
volume = new(max(t_var)+1, "integer")
start_time = new(max(t_var)+1, "integer")
end_time = new(max(t_var)+1, "integer")
delta_time = new(max(t_var)+1, "integer")
centroid_lat = new(max(t_var)+1, "double")
centroid_lon = new(max(t_var)+1, "double")
area = new(max(t_var)+1, "float", -9999.)
v_av = new(max(t_var)+1, "float", -9999.)
p_10 = new(max(t_var)+1, "float", -9999.)
p_25 = new(max(t_var)+1, "float", -9999.)
p_50 = new(max(t_var)+1, "float", -9999.)
p_75 = new(max(t_var)+1, "float", -9999.)
p_90 = new(max(t_var)+1, "float", -9999.)
t1D = ndtooned(t_var)
p1D = ndtooned(p_var)
dsizes_t = dimsizes(t_var)
do i=1,max(t_var)
indices_t = ind_resolve(ind(t1D.eq.i),dsizes_t)
volume(i) = num(t_var.eq.i)
start_time(i) = indices_t(0,0)
end_time(i) = indices_t(dimsizes(indices_t(:,0))-1,0)
delta_time(i) = 1+end_time(i)-start_time(i)
centroid_lat(i) = avg(lat(indices_t(:,1)))
centroid_lon(i) = avg(lon(indices_t(:,2)))
area(i) = volume(i)/delta_time(i)
v_av(i) = avg(p1D(ind(t1D.eq.i)))
p_10(i) = Percentile(p1D(ind(t1D.eq.i)),10)
p_25(i) = Percentile(p1D(ind(t1D.eq.i)),25)
p_50(i) = Percentile(p1D(ind(t1D.eq.i)),50)
p_75(i) = Percentile(p1D(ind(t1D.eq.i)),75)
p_90(i) = Percentile(p1D(ind(t1D.eq.i)),90)
delete(indices_t)
end do
; Write data as table to text file
r = ispan(1,max(t_var),1)
system("/bin/rm -f var.txt")
fname = "var.txt"
fhead = systemfunc("echo -e tnum $'\t' start $'\t' end $'\t' dt $'\t' c_lat $'\t' c_lon $'\t' vol $'\t' area $'\t' v_avg $'\t' p_10 $'\t' p_25 $'\t' p_50 $'\t' p_75 $'\t' p_90 >> "+fname)
print(fhead)
do i=1,max(t_var)
str_var = sprinti("%8.0i",r(i-1))+"$'\t'"+sprinti("%4.0i",start_time(i))+"$'\t'"+sprinti("%4.0i",end_time(i))+"$'\t'"+sprinti("%4.0i",delta_time(i))+"$'\t'"+\
sprintf("%2.2f",centroid_lat(i))+"$'\t'"+sprintf("%3.2f",centroid_lon(i))+"$'\t'"+\
sprinti("%10.0i",volume(i))+"$'\t'"+sprintf("%8.2f",area(i))+"$'\t'"+sprintf("%3.2f",v_av(i))+"$'\t'"+\
sprintf("%3.2f",p_10(i))+"$'\t'"+sprintf("%3.2f",p_25(i))+"$'\t'"+\
sprintf("%3.2f",p_50(i))+"$'\t'"+sprintf("%3.2f",p_75(i))+"$'\t'"+\
sprintf("%3.2f",p_90(i))
cmd = systemfunc("echo -e " + str_var + " >> "+fname)
print(cmd)
end do
print("Total run time: " + (get_cpu_time() - begTime))
end