NCL从全球格点数据提取中国区域数据的疑问

LeoBao

;*****************************************************
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 "./zh.ncl"
;****a****************************************************
begin

    f_BCC=addfile("/mnt/xfs/CMIP5/output/CMIP5/pr/historical/BCC-CSM/pr_Amon_BCC-CSM_historical_r1i1p1_185001-200512.nc","r")

  time=f_BCC->time
   lat=f_BCC->lat({5:65})
   lon=f_BCC->lon({65:145})
   nlat=dimsizes(lat)
   nlon=dimsizes(lon)
   ntime=dimsizes(time)
   
   FD=f_BCC->tas
   
;----------------the China area-------------------------

nx=701
ny=401
sx=70.0
ex=140.0
sy=15.0
ey=55.0

     stn=nlat*nlon
     lat2d=new((/nlat,nlon/),"float")
     lon2d=new((/nlat,nlon/),"float")

  do i=0,nlon-1
      lat2d(:,i)=dble2flt(lat)
    end do
  do j=0,nlat-1
     lon2d(j,:)=dble2flt(lon)
   end do
    lat1d=ndtooned(lat2d)
    lon1d=ndtooned(lon2d)
    x=new(5000,"float")
    y=new(5000,"float")
    tn=new((/5000/),"float")

    FD_CN=new((/45,401,701/),"float")
    delete(lat)
    delete(lon)

   lat=ispan(150,550,1)*0.1
   lon=ispan(700,1400,1)*0.1

    year!0="year"
    FD_CN!0="year"
    FD_CN!1="lat"
    FD_CN!2="lon"

    FD_CN&lat=lat
    FD_CN&lon=lon
    FD_CN&lat@units= "degrees_north"
    FD_CN&lon@units= "degrees_east"
   
    ts_reorde=FD(lat|:,lon|:,year|:)

    do t=0,44
    j=0
   nn=0
   tt=ndtooned(ts_reorde(:,:,t))
  do i=0,stn-1
  if(.not.(ismissing(tt(i))))then
   x(j)=lon1d(i)
   y(j)=lat1d(i)
   tn(j)=tt(i)
   j=j+1
  end if
end do
j=j-1
  FD_CN(t,:,:)=zh_rectgrid(x(0:j),y(0:j),tn(0:j),1.0,"./cnmask.dat","float",nx,ny,sx,ex,sy,ey,"nat",True)
  end do

  FD_T=new((/45/),float)
  pi=3.14159
  wgt=cos(lat/180*pi)
  FD_T=dim_avg_n(dim_avg_wgt_n(FD_CN,wgt,1,1),1)
  print(max(FD_T))
  FD_avg=avg(FD_T)
  print(FD_avg)

;------write the results-----------------------
opt = True
opt@fout = "FD_BCC.txt"
write_matrix(onedtond(FD_T,(/45,1/)),"f12.5", opt)
trendobs=regline(year,FD_T)
  print(trendobs)   
   end

zh.ncl中一段上面程序用到的定义的函数
function zh_rectgrid(tx:numeric, ty:numeric, tz:numeric, mvalue:numeric, mfile:string, mtype:string, nx:integer, ny:integer, sx:numeric, ex:numeric, sy:numeric, ey:numeric, mode:string, dataintaiwan:logical)
local gx,gy,cnmask,gdata,mdata,is,ie,js,je,opt
begin
  ; in cnmask.dat ( China: 1, Out: -1 )
        ; nx=701
        ; ny=401
        ; sx=70.0
        ; ex=140.0
        ; sy=15.0
        ; ey=55.0
  gx = fspan(sx, ex, nx)
  gx!0 = "lon"
  gx@units = "degrees_east"
  gx&lon = gx

  gy = fspan(sy, ey, ny)
  gy!0 = "lat"
  gy@units = "degrees_north"
  gy&lat = gy

  setfileoption("bin","ReadByteOrder","LittleEndian")
  cnmask = fbindirread(mfile, 0, (/ny,nx/), mtype)
  cnmask!0 = "lat"
  cnmask&lat = gy
  cnmask!1 = "lon"
  cnmask&lon = gx
  
  if (mode .eq. "nat") then
    gdata = natgrids(ty, tx, tz, gy, gx)
  end if
  if (mode .eq. "tri") then
    opt=True
    opt@distmx=100.0
    gdata = triple2grid(tx, ty, tz, gx, gy, opt)
  end if
  if (mode .eq. "css") then
    gdata = cssgrid_Wrap(ty, tx, tz, gy, gx)
  end if
  if (mode .eq. "csa") then
    knots = (/4,4/)
    gdata = csa2(ty, tx, tz, knots, gy, gx)
  end if
  if (mode .eq. "ds") then
    dssetp("dmx", 1.0)
    gdata = dsgrid2(ty, tx, tz, gy, gx)
  end if
  gdata!0 = "lat"
  gdata&lat = gy
  gdata!1 = "lon"
  gdata&lon = gx

  mdata = mask(gdata, cnmask, mvalue)
  mdata!0 = "lat"
  mdata&lat = gy
  mdata!1 = "lon"
  mdata&lon = gx  
  
  ; if there is not data in taiwan
  if (.not. dataintaiwan) then
    is=floattoint((120.0-sx)/0.1)
    ie=floattoint((122.0-sx)/0.1)
    js=floattoint((22.0-sy)/0.1)
    je=floattoint((25.1-sy)/0.1)
    mdata(js:je,is:ie)=mdata@_FillValue
  end if
  return mdata
end


这是我找到的提取中国区域的程序，但有一部分不太明白，
stn=nlat*nlon这个值是不是赋错了啊，不然下面的循环判断是不是缺测值的那段只执行了一部分点，因为做了reorder
ts_reorde=FD(lat|:,lon|:,year|:)

    do t=0,44
    j=0
   nn=0
   tt=ndtooned(ts_reorde(:,:,t))
  do i=0,stn-1
  if(.not.(ismissing(tt(i))))then
   x(j)=lon1d(i)
   y(j)=lat1d(i)
   tn(j)=tt(i)
   j=j+1
  end if
end do
j=j-1

矮油喂

vay 发表于 2019-5-29 21:36
有没有课代表解释下啥意思！我有好多全球nc数据文件不知道怎么提取中国的

请问你解决了吗全球格点数据怎么提取中国区域的，我也不会求解答

宋秀鹏

什么数据？

流浪的心

只有脚本，没有数据文件！

流浪的心

提供数据文件更好！

zhaoshy

楼主，请问你那里有cnmaster.dat文件吗？或者知道在哪里可以找到这个文件吗？万分感谢！

983019071

哪里有cnmaster.dat文件吗？    急

流浪的心

同求，求大神指点！

流浪的心

求答案！！！

vay

有没有课代表解释下啥意思！我有好多全球nc数据文件不知道怎么提取中国的