批量读取逐日数据并整合

threetee · 发表于 2014-9-24 11:55:12

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

有文件名如图的一些逐日数据，从1992年到2012年的，想用ncl批量读取然后整合成一个三维的数据，时间做一维。
请问下要怎么弄？谢谢

longlivehj · 发表于 2014-9-24 13:32:37

全部整起来有近30GB！

threetee · 发表于 2014-9-24 13:54:57

longlivehj 发表于 2014-9-24 13:32
全部整起来有近30GB！

我是想把这些数据处理成月平均的数据，之前看别人的帖子推荐了一个函数。我之前用matlab的程序处理了一下，感觉不太对，而且最后再ncl这边画图也画不出来，所以想试一下直接用ncl来处理看是什么效果。

言深深 · 发表于 2014-9-24 14:51:28

我处理过26G的，由于我不会自己读取nc文件，所以用grads转换成dat文件，然后fortran处理的。鉴于一次性读取这么大数据速度比较慢，所以按照时间，层处理成很多小文件。

longlivehj · 发表于 2014-9-24 15:53:54

threetee 发表于 2014-9-24 13:54
我是想把这些数据处理成月平均的数据，之前看别人的帖子推荐了一个函数。我之前用matlab的程序处理了一下 ...

逐年处理成月平均，然后再合并吧，这样数据量会小很多！
方法有多种，ncl可以（addfiles，calculate_monthly_values），cdo也可以（cdo cat，cdo monmean）。

threetee · 发表于 2014-9-24 17:04:52

longlivehj 发表于 2014-9-24 15:53
逐年处理成月平均，然后再合并吧，这样数据量会小很多！
方法有多种，ncl可以（addfiles，calculate_mon ...

我去看了这个函数addfiles，它给出的例子好像都是数字变化类似于从1到100这样子变得，我这个年月日都在变，几个例子都试了一遍，都没有办法读进去。我又回去用matlab重新处理了一下，还是一样画不出图来，而且我用这个数据算的流函数画图，在ncl里边定义的颜色好像没起作用。不知道到底是数据的问题还是程序得问题了。

longlivehj · 发表于 2014-9-24 17:16:44

threetee 发表于 2014-9-24 17:04
我去看了这个函数addfiles，它给出的例子好像都是数字变化类似于从1到100这样子变得，我这个年月日都在变 ...

不明白你讲的啥意思。
我既没有你的数据，也没有你的程序，所以不知道“到底是数据问题还是程序问题”。

threetee · 发表于 2014-9-24 20:42:36

longlivehj 发表于 2014-9-24 17:16
不明白你讲的啥意思。
我既没有你的数据，也没有你的程序，所以不知道“到底是数据问题还是程序问题”。

ecco_wind_curl.nc (26.66 MB, 下载次数: 0)

2014-9-24 20:38 上传

点击文件名下载附件
这是我的数据，程序是这样的。麻烦帮看一下
;********************************************************
;  Load NCL scripts
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl"

begin
;files = systemfunc("ls ./windstress*.nc") ;批量处理
; f    =addfile("windstress01.nc","r")    ;以只读方式打开一个文件
  f0 =addfile("./ecco_wind_curl.nc","r")
print(f0)

taux1    = f0->taux_mean(time|:,lat|:,lon|:)
tauy1    = f0->tauy_mean(time|:,lat|:,lon|:)
rota    = f0->rota_mean(time|:,lat|:,lon|:)
rota=rota*1e3
if (any(isnan_ieee(rota))) then       ;isnan_ieee返回包含ieee nan的数组
   value  = 1.e20
   replace_ieeenan (rota, value, 0) ;该函数用于将nan值转换成用户指定的值 ;value是用于代替nan
   rota@_FillValue = value
end if
if (any(isnan_ieee(taux1))) then       ;isnan_ieee返回包含ieee nan的数组
   value  = 1.e20
   replace_ieeenan (taux1, value, 0) ;该函数用于将nan值转换成用户指定的值 ;value是用于代替nan
   taux1@_FillValue = value
end if
if (any(isnan_ieee(tauy1))) then       ;isnan_ieee返回包含ieee nan的数组
   value  = 1.e20
   replace_ieeenan (tauy1, value, 0) ;该函数用于将nan值转换成用户指定的值 ;value是用于代替nan
   tauy1@_FillValue = value
end if

printVarSummary(taux1)

lat    =f0->lat
lon    =f0->lon
lat@units = "degrees_north"
lon@units = "degrees_east"

taux1_mean=dim_avg_n_Wrap(taux1,0)    ;求年平均
tauy1_mean=dim_avg_n_Wrap(tauy1,0)
rota_mean =dim_avg_n_Wrap(rota,0)
taux2    =dim_rmvmean_n_Wrap(taux1,0) ;直接是三维场时减去年平均
tauy2    =dim_rmvmean_n_Wrap(tauy1,0)
rota1    =dim_rmvmean_n_Wrap(rota,0)
printVarSummary(taux2)
taux1&lat@units = "degrees_north"
taux1&lon@units = "degrees_east"
tauy1&lat@units = "degrees_north"
tauy1&lon@units = "degrees_east"
rota&lat@units = "degrees_north"
rota&lon@units = "degrees_east"
; taux2&lat@units = "degrees_north"
; taux2&lon@units = "degrees_east"
; tauy2&lat@units = "degrees_north"
; tauy2&lon@units = "degrees_east"
rota1&lat@units = "degrees_north"
rota1&lon@units = "degrees_east"

taux2&lat=lat
taux2&lon=lon

  printVarSummary(taux2)
  delete(taux2@long_name)
  delete(tauy2@long_name)
  delete(rota1@long_name)
;printVarSummary(lat)
;******************************************
;---Graphics
  wks1             = gsn_open_wks("ps","ecco_wind_field_mon_mean_minus_annual_mean")
;  wks2             = gsn_open_wks("ps","ecco_wind_field_annual_mean111")

  gsn_define_colormap(wks1,"BlWhRe")
;  gsn_define_colormap(wks2,"BlWhRe")

  vres                = True
  plot             = new(12,graphic)

  vres@gsnDraw       = False
  vres@gsnFrame       = False
  vres@gsnAddCyclic = False ;只有数据是全球的才需要设置成True

  vres@vfXArray       = lon
  vres@vfYArray       = lat
  vres@sfXArray       = lon
  vres@sfYArray       = lat

  vres@cnFillOn          = True             ; turn on color for contours
  vres@cnLinesOn          = False             ; turn off contour lines
  vres@cnLineLabelsOn    = False
  vres@gsnScalarContour    = True
  vres@gsnSpreadColors    = True
  vres@mpMinLonF          = min(lon)
  vres@mpMaxLonF          = max(lon)
  vres@mpMinLatF          = min(lat)
  vres@mpMaxLatF          = max(lat)
  vres@lbLabelBarOn       = False

;*******设置contour的属性**********************************************************

vres@cnLevelSelectionMode ="ExplicitLevels"  ;设置colorbar等距
; vres@cnMinLevelValF       =-10
; vres@cnMaxLevelValF       =10
; vres@cnLevelSpacingF       =1
vres@cnLevels             =(/-10,-8,-6,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.7,-0.4,-0.2,-0.1,0,0.1,0.2,0.4,0.7,1,1.5,2,2.5,3,3.5,4,6,8,10/)

;*******设置vector的属性***********************************************************
  vres@vcRefMagnitudeF       = 0.03 ;设置矢量标度
  vres@vcRefLengthF          = 0.02
  vres@vcMinDistanceF          =0.04    ; 设置向量的疏密
  vres@vcRefAnnoOrthogonalPosF = -1.0
  vres@vcLineArrowThicknessF  =0.7
  vres@vcFillArrowWidthF    =0.03
  vres@vcVectorDrawOrder       = "PostDraw"

;******每个图形的大小***************************************************************
;  vres@vpHeightF                =0.4
;  vres@vpWidthF                =0.4

;******各图形坐标轴属性*************************************************************
  vres@tmYRLabelsOn       = False             ; no right labels
  vres@tmYROn             = False             ; no right tickmarks
  vres@tmXTLabelsOn       = False             ; do not draw top labels
  vres@tmXTOn             = False
  vres@tmXBOn             = False
  vres@tmXBLabelsOn       = False
  vres@tmXBOn             = False
  vres@tmYLLabelsOn       = False
  vres@tmYLOn             = False
;  print(1)
  do i=0,1
plot(i+1)=gsn_csm_vector_scalar_map_ce(wks1,taux2(i+1,:,:),tauy2(i+1,:,:),rota1(i+1,:,:),vres)
print(1)
plot(i+4)=gsn_csm_vector_scalar_map_ce(wks1,taux2(i+4,:,:),tauy2(i+4,:,:),rota1(i+4,:,:).vres)
print(1)
plot(i+7)=gsn_csm_vector_scalar_map_ce(wks1,taux2(i+7,:,:),tauy2(i+7,:,:),rota1(i+7,:,:),vres)
  end do
  print(1)
  vres@tmYLOn                =False
  vres@tmXBOn                =True
  vres@tmXBLabelsOn          =True
  vres@tmXBMode             ="Explicit"
  vres@tmXBValues          =(/25,40,55,70,85,100,115/)
  vres@tmXBLabels          =(/"25E","40E","55E","70E","85E","100E","115E"/)

  do i=0,1
plot(i+10)=gsn_csm_vector_scalar_map_ce(wks1,taux2(i+10,:,:),tauy2(i+10,:,:),rota1(i+10,:,:),vres)
  end do

  vres@tmYLOn                =True
  vres@tmYLLabelsOn          =True
  vres@tmYLMode             ="Explicit"
  vres@tmYLValues          =(/-30,-20,-10,0,10,20/)
  vres@tmYLLabels          =(/"30S","20S","10S","Eq","10N","20N"/)

  plot(9) =gsn_csm_vector_scalar_map_ce(wks1,taux2(9,:,:),tauy2(9,:,:),rota1(9,:,:),vres)

  vres@tmXBOn             =False

  plot(0) =gsn_csm_vector_scalar_map_ce(wks1,taux2(0,:,:),tauy2(0,:,:),rota1(0,:,:),vres)
  plot(3) =gsn_csm_vector_scalar_map_ce(wks1,taux2(3,:,:),tauy2(3,:,:),rota1(3,:,:),vres)
  plot(6) =gsn_csm_vector_scalar_map_ce(wks1,taux2(6,:,:),tauy2(6,:,:),rota1(6,:,:),vres)

;panel plot ;所有图画在一张图上
pres             =True

pres@gsnPanelFigureStrings= (/"JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG","SEP","OCT","NOV","DEC"/)
pres@amJust = "BottomLeft"
; pres@gsnPanelScalePlotIndex       = 3
; pres@pmLabelBarOrthogonalPosF       =0.005
pres@gsnPanelFigureStringsFontHeightF=0.01
pres@gsnPanelXWhiteSpacePercent    =0.0
pres@gsnPanelYWhiteSpacePercent    =0
pres@gsnPanelLabelBar       =True
pres@lbOrientation    ="Vertical"
pres@pmLabelBarHeightF    =0.5
pres@lbLabelStride       = 1
pres@gsnPanelScalePlotIndex    = 7
pres@gsnPanelLeft             = 0.1
pres@gsnPanelRight             = 0.95
;  drawNDCGrid(wks1)
pres@txFontHeightF = 0.01
gsn_text_ndc(wks1,"x10~S~-3",0.92,0.76,pres)

gsn_panel(wks1,plot,(/4,3/),pres)

vres@tmXBOn       =True
vres@lbLabelBarOn =True
plot2             =gsn_csm_vector_scalar_map_ce(wks2,taux1_mean,tauy1_mean,rota_mean,vres)
  end

fzghui · 发表于 2014-9-24 21:04:38

。。。。好复杂

longlivehj · 发表于 2014-9-24 21:40:41

threetee 发表于 2014-9-24 20:42
这是我的数据，程序是这样的。麻烦帮看一下
;***************************************************** ...

跪了！
说好的“批量读取逐日数据并整合”呢！

		自动登录	找回密码
密码			立即注册

[其他] 批量读取逐日数据并整合

登录后查看更多精彩内容~