中国气象科学数据0.5°×0.5°格点数据的处理小结

徵言片羽

    科研小白，白嫖多年，第一次发帖，经验不足，还请各路大神指正。话不多说，直入正题。
    最近看论文，想重现一下文中降水的年际变化，需要用到中国气象科学数据网所提供0.5°×0.5°格点降水月数据。下载过程非常曲折，首先要实名注册，通过后才能下载，并且等级较低时每天只能下载20M数据，我签到了好多天，才能下得稍微多一点。（网上不少人吐槽这个事，咱自己用数据还这么费劲）
   格点数据是文本格式的，并且是单独的月数据，我下载了51年数据，共612个txt数据，真是崩溃！在txt中打开以后，如果缩小后，会发现数据是个中国地图的轮廓。但不含有经纬度信息，ncl不能根据经纬度信息，划定所需地区的范围，也就无法直接拿来画图，所以我需要把612个数据转换为netCDF数据，也就是常用的.nc数据，一共三个维度：time,lat,lon，lat和lon可以根据数据格式说明，来反向还原，其中关键的语句是：
ncols                = 128
nrows                = 72
xllcorner        = 72.25  ;lon
yllcorner        = 18.25  ;lat
lat                  = fspan (yllcorner,yllcorner+(nrows-1)*0.5,nrows)
lon                   = fspan (xllcorner,xllcorner+(ncols-1)*0.5,ncols)
lat@units      = "degrees_north"
lon@units      = "degrees-east"

时间维度还没解决，等解决了来完善此帖。完整代码如下：
begin
  diri                = "D:/data/data_create/pre-nc/"
  filname         = "SURF_CLI_CHN_PRE_MON_GRID_0.5-"
  filformat        = ".nc"

  files        = systemfunc("ls /cygdrive/d/data/pre/SURF_*")            ; ?????????
  nf                = dimsizes(files)
  do i=1961,2011
   do j=1,12
    k=(i-1961)*12+j-1                                                        ;??grid?????
    print(i*100+j +" Now reading... " + files(k))     ;????
        filtime = i*100+j
        outfile         = diri+filname+filtime+filformat       
       
;************??nc???????*************************
system("rm -f "+outfile)
ncdf                 = addfile(outfile,"c")
dim_names         = (/ "lat",  "lon" /)
dim_Unlim         = (/ False , False/)
dim_sizes         = (/ 72,128/)
filedimdef(ncdf,dim_names,dim_sizes,dim_Unlim)
ncols                = 128
nrows                = 72
xllcorner        = 72.25  ;lon
yllcorner        = 18.25  ;lat
lat                  = fspan (yllcorner,yllcorner+(nrows-1)*0.5,nrows)
lon                   = fspan (xllcorner,xllcorner+(ncols-1)*0.5,ncols)
lat@units      = "degrees_north"
lon@units      = "degrees-east"
;********************?grid????????nc???********************
T = readAsciiTable(files(k),ncols,"float",6)
p= T(::-1,:)                       
p@_FillValue        = -9999.0
p!0            = "lat"                     
p!1                   = "lon"
p&lat          = lat                    
p&lon          = lon
p@long_name        = "Monthly mean precipition in China"
p@units            = "mm"
ncdf->p          = p

delete(filtime)               ;?????????????
delete(outfile)
delete(T)
delete(p)
delete(ncdf)

end do                        
end do                              


end

气象小鬼

这个降水数据在有的时间里都没数据