用NCL转换nc数据格式

zzzzzx

这是nc数据的头文件：
netcdf mh_results_18_27.daymax {
dimensions:
    x = 163 ;
    y = 127 ;
    time = UNLIMITED ; // (3652 currently)
variables:
    double time(time) ;
        time:standard_name = "time" ;
        time:long_name = "Time" ;
        time:units = "hours since 1901-01-01 00:00:00" ;
        time:calendar = "standard" ;
    float T2(time, y, x) ;
        T2:units = "K" ;
        T2:FieldType = 104 ;
        T2:MemoryOrder = "XY " ;
        T2:description = "TEMP at 2 M" ;
    float TSK(time, y, x) ;
        TSK:units = "K" ;
        TSK:FieldType = 104 ;
        TSK:MemoryOrder = "XY " ;
        TSK:description = "SURFACE SKIN TEMPERATURE" ;

// global attributes:
        :CDI = "Climate Data Interface version 1.6.0 (http://code.zmaw.de/projects/cdi)" ;
        :Conventions = "CF-1.4" ;
        :history = "Tue Jul 01 12:01:34 2014: cdo daymax mh_results_18_27_TSK.nc mh_results_18_27.daymax.nc\n",
            "Tue Jul 01 11:57:40 2014: cdo selname,T2,TSK mh_results_18_27.nc mh_results_18_27_TSK.nc\n",
            "Tue Jul 01 11:52:03 2014: cdo selyear,18/27 mh_results.nc mh_results_18_27.nc\n",
            "Tue Jul 01 11:43:29 2014: cdo copy wrf_mh_t2_0016.nc wrf_mh_t2_0017.nc wrf_mh_t2_0018.nc wrf_mh_t2_0019.nc wrf_mh_t2_0020.nc wrf_mh_t2_0021.nc wrf_mh_t2_0022.nc wrf_mh_t2_0023.nc wrf_mh_t2_0024.nc wrf_mh_t2_0025.nc wrf_mh_t2_0026.nc wrf_mh_t2_0027.nc wrf_mh_t2_0028.nc mh_results.nc" ;
        :CDO = "Climate Data Operators version 1.6.0 (http://code.zmaw.de/projects/cdo)" ;
}

longlivehj

zzzzzx 发表于 2014-7-15 20:39
找人要了个相同区域的文件，这次有了...
        float xlat(y, x) ;
                xlat:units = "d ...

DefineMode使用了中文的引号！

longlivehj

转换成什么格式？
没有lat、lon，只有x、y，也得想办法解释啊！通过空间范围、分辨率等等推算。

zzzzzx

longlivehj 发表于 2014-7-14 22:02
转换成什么格式？
没有lat、lon，只有x、y，也得想办法解释啊！通过空间范围、分辨率等等推算。

我是想转换成符合COARDS格式的nc文件，或者是dat文件然后自己编ctl。
只是分辨率不知道，从头文件里只能看到x y的格点数，起始位置也不知道，这改怎么做呢

longlivehj

zzzzzx 发表于 2014-7-15 13:22
我是想转换成符合COARDS格式的nc文件，或者是dat文件然后自己编ctl。
只是分辨率不知道，从头文件里只能 ...

数据从什么地方来的？
从global attributes中的history来看，数据貌似跟wrf有关，用cdo进行过一系列处理。这些处理是不是你自己做的？

zzzzzx

longlivehj 发表于 2014-7-15 14:20
数据从什么地方来的？
从global attributes中的history来看，数据貌似跟wrf有关，用cdo进行过一系列处理 ...

貌似是wrf模拟出来的数据，全国范围的。不是我处理的，不然wrfout就不会有那么麻烦了...甚至这都不是我用的，我帮着处理一下转换成grads能用的格式给别人画图...
也就是说光凭这个文件的头文件信息确实不能确定经纬度信息了吧？要是这样我再去问问这资料来源

longlivehj

zzzzzx 发表于 2014-7-15 14:26
貌似是wrf模拟出来的数据，全国范围的。不是我处理的，不然wrfout就不会有那么麻烦了...甚至这都不是我用 ...

嗯，是的，netcdf本来是自解释的文件，元数据充足。但是，数据处理过程中，难免有信息丢失。只能从源头找起了。

zzzzzx

longlivehj 发表于 2014-7-15 14:46
嗯，是的，netcdf本来是自解释的文件，元数据充足。但是，数据处理过程中，难免有信息丢失。只能从源头找 ...

恩，好的。很感谢你的耐心解答！

zzzzzx

longlivehj 发表于 2014-7-15 14:46
嗯，是的，netcdf本来是自解释的文件，元数据充足。但是，数据处理过程中，难免有信息丢失。只能从源头找 ...

找人要了个相同区域的文件，这次有了...
        float xlat(y, x) ;
                xlat:units = "degree_north" ;
                xlat:FieldType = 104 ;
                xlat:MemoryOrder = "XY " ;
                xlat:description = "LATITUDE, SOUTH IS NEGATIVE" ;
        float xlong(y, x) ;
                xlong:units = "degree_east" ;
                xlong:FieldType = 104 ;
                xlong:MemoryOrder = "XY " ;
                xlong:description = "LONGITUDE, WEST IS NEGATIVE" ;
然后我想创建个新的nc文件，然后把time,lat,lon,变量根据原文件直接赋值，这样让新创建的nc文件能直接被grads打开，这样可行吗？我试了下，东拼西凑的ncl脚本有问题...

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
begin
;---read in a netCDF variable
   fi   = addfile ("pi_results_18_27.daymean.nc", "r")
   lat  = fi->xlat                        ; read in netCDF variable
   time = fi->time
   lon  = fi->xlong
   t2   = fi->T2

   fout = addfile("out.nc","c")
   
   ; create global attributes
   fileAtt                = True
   fileAtt@title          = "Sample"
   fileAtt@Conventions    = "None"
   fileAtt@creation_date  = systemfunc("date")
   setfileoption(fout,“DefineMode”,True)        ; optional
   fileattdef(fout,fileAtt)

   ; predefine coordinate variables
   dimNames = (/"time","lat","lon"/)
   dimSizes = (/-1,nlat,nlon/)   ; -1 means unspecified
   dimUnlim = (/True,False,False/)

   ; predefine names, type, dimensions
   ; explicit dimension naming or getvardims can be used
   filedimdef(fout,dimNames,dimSizes,dimUnlim)
   filevardef(fout,"time",typeof(time),getvardims(time))
   filevardef(fout,"lat" ,typeof(lat) ,"lat")
   filevardef(fout,"lon" ,"float" ,"lon")
   filevardef(fout,"TMP" ,typeof(T) , getvardims(T) )

   ; predefine each variable’s attributes
   filevarattdef(fout,"time",time)
   filevarattdef(fout,"lat" ,lat)
   filevarattdef(fout,"lon" ,lon)
   filevarattdef(fout,"t2" ,t2)
   setfileoption(fout,”SuppressDefineMode”,False) ; optional

   ; output values only
   fout->time    =  (/time/)
   fout->lat     =  (/lat/)
   fout->lon     =  (/lon/)
   fout->t2      =  (/t2/)  ;t2

end

初学ncl,问题估计很多...能否请你帮忙看一下呢?错误提示是
fatal:syntax error: line 18 in file gaigeshi.ncl before or near setfileoption(fout,

zzzzzx

longlivehj 发表于 2014-7-15 14:46
嗯，是的，netcdf本来是自解释的文件，元数据充足。但是，数据处理过程中，难免有信息丢失。只能从源头找 ...

顺便再问个比较小白的问题...像这种地理坐标是2维的,
lat=fi->xlat
lon=fi->xlong  
data=fi->t2
data@lat2d=lat
data@lon2d=lon
用上面几句可以给数据附加经纬度信息吗？那我输出变量时经纬度信息也会输出吗？