首页 > 解决方案 > 是否有 python 函数来获取 netcdf 索引、按值计数等

问题描述

我在 python 方面没有经验,我正在将 NCL 脚本转换为 python,希望 python 运行得更快。环顾四周,我没有找到我认为是 NCL 脚本中最简单计算的答案。看看如何完成更艰难的计算,我也没有找到如何在 python 中完成这些的答案。

大部分计算是在将 3 维变量转换为 1 维变量并查询它们在数组空间中的值和位置之后完成的。知道了 t 个变量在数组空间中的位置,我们可以得到对应于 t 个变量整数值的 p 个变量值。

计算如下:

所有这些计算都将值保存在一维数组中,维度大小等于 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

标签: pythonnetcdfncl

解决方案


可能存储在 NetCDF 文件中的数据种类繁多,我无法向您保证它会解决您的问题,但是如果您需要,GDAL库的 Python 绑定及其专用的 NetCDF 驱动程序肯定是一个很好的起点通过 Python 代码探索这种格式的文件。

一旦您能够使用此库访问您的数据,如果您准备将其分解成更小的部分而不是将整个 NCL 脚本扔在那里进行转换,您很可能会在这里找到解决实际问题的帮助进入 Python。


推荐阅读