NCL处理数据

四月雨之音

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"
begin

;----------------------------------------  GPCP -----------------------------
dlonL = 60
dlonR = 160
dlatS = -30
dlatN = 30
nlon=100
nlat=60

f1 = addfile("/home/yangsong3/data-observation/PRECPT/gpcp_1dd_v1.2_p1d.19961001-20150531.daily.nc","r")
tc_date = ut_calendar(f1&time,-2)
print(tc_date)
print(tc_date(0:1))
pStart  = 20000101
pEnd    = 20141231
tstart  = ind(tc_date.eq.pStart)
tend    = ind(tc_date.eq.pEnd)

temp = f1->precip(tstart:tend,::-1,:)
gpcp = temp(0:5474,{dlatS:dlatN}, {dlonL:dlonR})

mmm=0
do tt = 0,5478
if(tt.ne.59.and.tt.ne.1520.and.tt.ne.2981.and.tt.ne.4442)
gpcp(mmm,{dlatS:dlatN}, {dlonL:dlonR})= temp(tt,{dlatS:dlatN}, {dlonL:dlonR})
mmm=mmm+1
end if
end do
print(mmm+"=5474?")
print(mmm)
year = 15   ; 2000-2014
day  = 90   ;Jan(31)+Feb(28)+Mar(31)
gpcp_JFM = new((/year,day,nlat,nlon/),"float")
do i=0,year-1
  do j=0,day-1
   gpcp_JFM(i,j,:,:)=gpcp(i*365+j,:,:)
   end do
end do
gpcp_JFM_mean =  dim_avg_n_Wrap(gpcp_JFM,1)

printVarSummary(gpcp_JFM)
printVarSummary(gpcp_JFM_mean)

;_____________________________




nlat = 190
nlon = 384

fall = new((/5475,nlat,nlon/),"float")

mmm = 0
do year = 2000,2014
all_files = systemfunc("ls /home/yangsong3/data-model/CFSv2_2000_2014/prate_daily/cfsv2.daily.pratesfc."+year+"*.01.dat")
   if(year.eq.2000.or.year.eq.2004.or.year.eq.2008.or.year.eq.2012)
      do ii = 0,dimsizes(all_files)-1
           if(ii.ne.59)
               fall(mmm,:,:) = fbindirread(all_files(ii),0,(/nlat,nlon/),"float")
               mmm = mmm+1
          end if
      end do  
   else
      do ii = 0,dimsizes(all_files)-1
         fall(mmm,:,:) = fbindirread(all_files(ii),0,(/nlat,nlon/),"float")
         mmm = mmm+1
     end do
  end if
  delete(all_files)
end do

fall = fall*3600*24
print("ok")
;_____interpolate CFSV to GPCP____________
nnlat = 180
nnlon = 360
cfsv2_pre  = new((/5475,nnlat,nnlon/),"float")
nlat_old  = (/-89.277, -88.340, -87.397, -86.454, -85.509, -84.565, -83.620, -82.676, -81.731, -80.786, -79.841, -78.897, -77.952, -77.007,\
-76.062, -75.117, -74.173, -73.228, -72.283, -71.338, -70.393, -69.448, -68.503, -67.559, -66.614, -65.669, -64.724, -63.779,\
  -62.834, -61.889, -60.945, -60.000, -59.055, -58.110, -57.165, -56.220, -55.275, -54.330,-53.386, -52.441, -51.496, -50.551,\
   -49.606, -48.661, -47.716, -46.771, -45.827, -44.882, -43.937, -42.992, -42.047, -41.102, -40.157, -39.212, -38.268, -37.323,\
    -36.378, -35.433, -34.488, -33.543, -32.598, -31.653, -30.709, -29.764, -28.819, -27.874, -26.929 ,-25.984, -25.039, -24.094,\
     -23.150, -22.205, -21.260, -20.315, -19.370, -18.425, -17.480, -16.535, -15.590, -14.646, -13.701, -12.756, -11.811, -10.866,\
       -9.921,  -8.976,  -8.031,  -7.087,  -6.142,  -5.197,  -4.252,  -3.307,  -2.362,  -1.417,  -0.472,   0.472,   1.417,   2.362,\
          3.307,   4.252,   5.197,   6.142,   7.087,   8.031,   8.976,   9.921,  10.866,  11.811,  12.756,  13.701,  14.646,  15.590,\
            16.535,  17.480,  18.425,  19.370,  20.315,  21.260,  22.205,  23.150,  24.094,  25.039,  25.984,  26.929,  27.874,  28.819,\
              29.764,  30.709,  31.653,  32.598,  33.543,  34.488,  35.433,  36.378,  37.323,  38.268,  39.212,  40.157,  41.102,  42.047,\
                42.992,  43.937,  44.882,  45.827,  46.771,  47.716,  48.661,  49.606,  50.551,  51.496,  52.441,  53.386,  54.330,  55.275,\
                  56.220,  57.165,  58.110,  59.055,  60.000,  60.945,  61.889,  62.834,  63.779,  64.724,  65.669,  66.614,  67.559,  68.503,\
                    69.448,  70.393,  71.338,  72.283,  73.228,  74.173,  75.117,  76.062,  77.007,  77.952,  78.897,  79.841,  80.786,  81.731,\
                      82.676,  83.620,  84.565,  85.509,  86.454,  87.397,  88.340,  89.277/)
nlon_old = fspan(0,359.0625,nlon)
nlat_new = fspan(-89.5,89.5,nnlat)
nlon_new = fspan(0.5,359.5,nnlon)
cfsv2_pre = linint2 (nlon_old,nlat_old,fall,True,nlon_new,nlat_new,0)

cfsv2_pre!0 = "time"

cfsv2_pre!1 = "lat"
cfsv2_pre!2 = "lon"
cfsv2_pre&lon  = fspan(0.5,359.5,nnlon)
cfsv2_pre&lat = fspan(-89.5,89.5,nnlat)
cfsv2_pre&time= fspan(0,5474,5475)

print("ok")
;_______________________________________________
nlon=100
nlat=60
year = 15   ; 2000-2014
day  = 90   ;Jan(31)+Feb(28)+Mar(31)
cfsv2_pre_JFM=new((/year,day,nlat,nlon/),"float")
do i=0,year-1
  do j=0,day-1
   cfsv2_pre_JFM(i,j,:,:)=cfsv2_pre(i*365+j,{dlatS:dlatN}, {dlonL:dlonR})
   end do
end do

cfsv2_pre_JFM!0 = "year"

cfsv2_pre_JFM!1 = "day"
cfsv2_pre_JFM!2 = "lat"
cfsv2_pre_JFM!3 = "lon"
cfsv2_pre_JFM&lon  = fspan(60.5,159.5,nlon)
cfsv2_pre_JFM&lat = fspan(-29.5,29.5,nlat)
cfsv2_pre_JFM&year  = fspan(2000,2014,year)
cfsv2_pre_JFM&day = fspan(0,89,day)

copy_VarCoords(cfsv2_pre_JFM, gpcp_JFM)

cfsv2_pre_JFM_mean = dim_avg_n_Wrap(cfsv2_pre_JFM,1)
printVarSummary(cfsv2_pre_JFM_mean)



cfsv2_pre_JFM_abn  = dim_rmvmean_n_Wrap(cfsv2_pre_JFM,1)



pre_cfsv_obs=new((/year,day+1,nlat,nlon/),"float")
do i=0,day-1
pre_cfsv_obs(:,i,:,:)=cfsv2_pre_JFM(:,i,:,:)
end do
pre_cfsv_obs(:,day,:,:)=gpcp_JFM_mean;（：，：，：）也可以  (:,0,:,:)error

print(pre_cfsv_obs(2,90,4,3))
print(gpcp_JFM_mean(2,4,3))
end

四月雨之音

loftydongshi 发表于 2016-7-18 09:22
没大懂你的意思，第九十一天的数值？如果三维的数组付给四维数组时已经定义好单点某个时刻就不会出现这种 ...

第一段是读取GPCP20000101到20141231的住日数据进来，但是把闰年的29日踢掉，
第二段是读取384*190的模式的数据，同样时间段剔除掉闰年29日
第三段是插值，将GPCP的作为最终插值的规格即360*180
第四段是节选60.5~159.5LON   -29.5~29.5LAT ,将数据处理成15年90天的规格，我只要1到3月一共九十天的数据
我最后要的数据有，90天平均后的15年的模式数据，90天鞠萍的15年的模式数据，90天模式的的外加GPCP90天平均后的15年的数据，有点复杂说的

十字天阶

longlivehj 发表于 2014-9-5 12:00
“挨着变量看，数值都正确；放到模式中，数据有错误”，到底是正确还是错误啊？错误的详情？
程序看上去没 ...

亲，你好，我这里用了三个do while循环，end do时一直报错，是不是也是因为NCL效率降低，罢工了啊？一般情况下，是不是尽量不要超过一个循环嵌套啊？谢谢亲