NCL linint2_Wrap 插值结果为空文件

龙哥的江小白

用ncl对cmip6 pr资料进行插值，想把1.875×1.25 插值成1×1°，利用linint2_Wrap插值得到的结果是个空文件~~不知道是哪里的问题，请大神指教~~
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
f = addfile("/mnt/d/bilinear/pr_day_ACCESS-ESM1-5_historical_r1i1p1f1_gn_19500101-19991231.nc","r")
pr = f->pr(:,:,:)
lon = f->lon(:)      ;(0-358.125) 1.875
lat = f->lat(:)       ;(-90-90)  1.25
time = f->time(:)  ;18262d
pr@_FillValue = 1.000000020040877e+20
pr@units = "kg m-2 s-1"
newlat = fspan(-90,90,181)
newlon = fspan(0,359,360)
newlat@units = "degrees_north"
newlon@units = "degrees_east"
newpr=new((/18262,181,360/),typeof(pr))
newpr = linint2_Wrap(lon,lat,pr,True,newlon,newlat,0)
newpr!0 = "time"
newpr!1 = "lat"
newpr!2 = "lon"
newpr&lat = newlat
newpr&lon = newlon
setfileoption("nc","Format","LargeFile")
fout = addfile("/mnt/d/bilinear/pr_day_ACCESS19500101-19991231.nc","c")

end

龙哥的江小白

运行过程中未提示错误，实在想不出是哪的问题

尽头的尽头

你只是创建了这个文件，你还没把变量写进去呢，当然是空的。
参考这里的两个方法http://www.ncl.ucar.edu/Applications/write_netcdf.shtml

雷子

用cdo吧，一条命令解决

龙哥的江小白

尽头的尽头 发表于 2020-12-10 19:14
你只是创建了这个文件，你还没把变量写进去呢，当然是空的。
参考这里的两个方法http://www.ncl.ucar.edu/ ...

是这样啊，学习了，谢谢！

龙哥的江小白

雷子 发表于 2020-12-10 19:57
用cdo吧，一条命令解决

请问是啥命令啊 这么厉害的！

雷子

龙哥的江小白 发表于 2020-12-10 20:07
请问是啥命令啊 这么厉害的！

cdo remapbil, r360*180 input.nc output.nc
remapbil是双线性插值，可改成其他

龙哥的江小白

尽头的尽头 发表于 2020-12-10 19:14
你只是创建了这个文件，你还没把变量写进去呢，当然是空的。
参考这里的两个方法http://www.ncl.ucar.edu/ ...

抱歉，不知道怎么上传图片，就把代码直接打上来了，前面的数字表示行号：
1 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl":
2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

6 begin
7 f = addfile("/mnt/d/bilinear/pr_day.nc","r")
8 pr = f->pr(:,:,:)
9 lon = f->lon(:)   ;(0-358.125) 1.875
10 lat = f->lat(:)   ;(-90-90)  1.25
11 time = f->time(:)  ;18262d
12 pr@_FillValue = 1e+20
13 pr@units = "kg m-2 s-1"
14 newlat = fspan(-90,90,181)
15 newlon = fspan(0,359,360)
16 newlat@units = "degrees_north"
17 newlon@units = "degrees_east"
18 newpr=new((/18262,181,360/),typeof(pr),1e+20)
19 newpr = linint2_Wrap(lon,lat,pr,True,newlon,newlat,0)
20 newpr!0 = "time"
21 newpr!1 = "lat"
22 newpr!2 = "lon"
23 newpr&lat = newlat
24 newpr&lon = newlon
25 ntim  = dimsizes(time)
26 nlat  = dimsizes(newlat)
27 nlon  = dimsizes(newlon)
28 setfileoption("nc","Format","LargeFile")
29 fout = addfile("/mnt/d/bilinear/pr_day_ACCESS19500101-19991231.nc","c")
30 setfileoption(fout,"DefineMode",True)
31 dimNames = (/"time", "lat", "lon"/)  
32 dimSizes = (/ ntim , nlat,  nlon/)
33 dimUnlim = (/ True , False, False/)   
34 filedimdef(fout,dimNames,dimSizes,dimUnlim)
35 filevardef(fout, "time" ,typeof(time),getvardims(time))                           
36 filevardef(fout, "lat"  ,typeof(lat),getvardims(lat))                          
37 filevardef(fout, "lon"  ,typeof(lon),getvardims(lon))                          
38 filevardef(fout, "pr"   ,typeof(pr) ,getvardims(pr))   
39 filevarattdef(fout,"pr",newpr)                           
40 filevarattdef(fout,"time" ,time)                    
41 filevarattdef(fout,"lat"  ,newlat)                     
42 filevarattdef(fout,"lon"  ,newlon)                     
43 fout->time   = (/time/)     
44 fout->lat    = (/newlat/)
45 fout->lon    = (/newlon/)
46 fout->pr     = (/newpr/)

48 end
运行过程中报错：fatal:Could not creat (/mnt/d/bilinear/pr_dayACCESS19500101-19991231.nc)
                         warning:FileSetFileOption:invalid file or format
  fatal:["Execute.c":8637]:Execute:Error occurred at or near line 34 in file /mnt/d/ncl/bilinear1211.ncl
不能创建输出文件，文件无效或者格式无效？不知道是怎么回事。第34行我是参照官网上写的，不知道为啥也报错了。可以麻烦您帮我看下嘛，谢谢！

龙哥的江小白

雷子 发表于 2020-12-11 09:13
cdo remapbil, r360*180 input.nc output.nc
remapbil是双线性插值，可改成其他

好的，谢谢啦，我试试！

wangchense

我觉得你这个输出文件有点麻烦，34行之后你可以参考一下我的这个
system("rm -f MODIS-annual.nc")
u=addfile("MODIS-annual.nc","c")
u@title="od550aer"
u->od550aer=model
其中model就是插值完的变量