【已解决】NCL站点插值函数报错

muggle

在同学的帮忙下，我写了一个统计高温日数的脚本，但是运行的时候在站点插值函数处报错，说左右数组维数不等，可是同学他那边同样的写法却可以运行，所以我是哪里错了。。。麻烦各位指点一下
报错结果如下
fatal:Dimension size mismatch on subscript #1, left-hand and right-hand side dimensions do not match
fatal:["Execute.c":8578]:Execute: Error occurred at or near line 88 in file /cygdrive/f/192/ncl/sta_num_of_tmax.ncl


下面是我的脚本

;----------------------------------------------------------------------
; sta_num_0f_tmax.ncl
;高温日数的统计及计算

;
; These files are loaded by default in NCL V6.2.0 and newer
; 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"

begin

;  读取站点逐日资料
  filename = "/cygdrive/f/192/1/EOF/data/35summer153dayobs.txt"

; Calculate the number of columns.
  ncols = numAsciiCol(filename)

; Given the # of columns, we can use readAsciiTable to read this file.
  data = readAsciiTable(filename,ncols,"float",0)

  nrows = dimsizes(data(:,0))    ; calculate # of rows
  print("'" + filename + "' has " + nrows + " rows and " + ncols + \
       " columns of data.")      

data_daily = new((/35,191,153/),"float")
do i=1,35
   do j=1,191
     data_daily ((i-1),(j-1),:)= data(((i-1)*191*153+(j-1)*153+0):((i-1)*191*153+(j-1)*153+152),4)
   end do
end do


;   printVarSummary(data)

;  读取阈值资料

   data_yz   = asciiread("/cygdrive/f/192/1/yuzhi/r95t35-54909.txt",(/35,191/),"float")
   data_yz!0   = "time"
   data_yz!1  =  "sta"
;   print(data_yz(0,190))

;     读取站点经纬度信息
   temp  = asciiread("/cygdrive/f/192/1/EOF/lon_lat-54909.txt",(/191,3/),"float")
   lat   = temp(:,2)
   lon   = temp(:,1)

;  print(lat(0))
;  print(lon(1))

;  判断大于阈值的高温日数
   num_tmax = new ((/35,191/),"float")
   do i=0,34
     do j=0,190
       num_tmax (i,j) =  num (data_daily(i,j,:) .gt. data_yz(i,j) )
     end do
   end do

  print(num_tmax(34,190))


;   构造格点   
   olon  =new(29,"float")
   olat  =new(17,"float")
   data1 =new((/35,17,29/),"float")
  do i=0,28
    olon(i)=70+i*2.5
  end do
  do j=0,16
    olat(j)=15+j*2.5
  end do


  olon!0           ="lon"
  olon@long_name   = "lon"
  olon@units       = "degrees_east"
  olon&lon         = olon

  olat!0           = "lat"
  olat@long_name   = "lat"
  olat@units       = "degrees_north"
  olat&lat         = olat


  rscan=(/10,5,3/)

  do i=0,34
      data1(i,:,:)  =  obj_anal_ic_Wrap(lon,lat,num_tmax(i,:),olon,olat,rscan,False)
  end do


end

muggle

画图部分由于还没运行到，所以就没给出来（其实我就是想看看自己回自己的帖子，能不能拿钱）

男紫汗

一般人不会从头看到尾吧，哪一行么标记或者高亮一下

muggle

男紫汗 发表于 2017-11-18 22:05
一般人不会从头看到尾吧，哪一行么标记或者高亮一下

恩恩，高亮了那一句了

muggle

我后来改为这样
  temp2 = new((/35,17,29/),"float")
  temp2(1,:,:) = obj_anal_ic_Wrap(lon,lat,num_tmax(1,:),olon,olat,rscan,False)
  printVarSummary(temp2)
报错为：
(0)     check_for_y_lat_coord: Warning: Data either does not contain a valid latitude coordinate array or doesn't contain one at all.
(0)     A valid latitude coordinate array should have a 'units' attribute equal to one of the following values:
(0)         'degrees_north' 'degrees-north' 'degree_north' 'degrees north' 'degrees_N' 'Degrees_north' 'degree_N' 'degreeN' 'degreesN' 'deg north'
(0)     check_for_lon_coord: Warning: Data either does not contain a valid longitude coordinate array or doesn't contain one at all.
(0)     A valid longitude coordinate array should have a 'units' attribute equal to one of the following values:
(0)         'degrees_east' 'degrees-east' 'degree_east' 'degrees east' 'degrees_E' 'Degrees_east' 'degree_E' 'degreeE' 'degreesE' 'deg east'
warning:ContourPlotInitialize: no valid values in scalar field; ContourPlot not possible:[errno=1101]

是不是意味着我的经纬度出了错

muggle

找到原因了
;     读取站点经纬度信息
   temp  = asciiread("/cygdrive/f/192/1/EOF/lon_lat-54909.txt",(/191,3/),"float")
   lat   = temp(:,2)
   lon   = temp(:,1)
  lon@units      = "degrees_east"
   lat@units      = "degrees_north"
这里得加上对站点经纬度单位的赋值

rabin_xu

用ncl来做这些计算的事情不是很慢吗？

muggle

rabin_xu 发表于 2017-11-20 12:49
用ncl来做这些计算的事情不是很慢吗？

是不如fortran快，但是我就是垂涎ncl的内置函数，调用起来太方便了。。

rabin_xu

muggle 发表于 2017-11-20 12:52
是不如fortran快，但是我就是垂涎ncl的内置函数，调用起来太方便了。。

曾经用hncl编了一个数据处理程序，具体忘了是算啥了，反正是循环计算的，数据量也比较大，是1040×680的吧，算的超慢，一算就是一天。后来只好自己又写fortran程序，眨眼就算完了。

muggle

rabin_xu 发表于 2017-11-20 12:56
曾经用hncl编了一个数据处理程序，具体忘了是算啥了，反正是循环计算的，数据量也比较大，是1040×680的 ...

你这个确实太大了。。目前我做的还比较简单，先用着试试，实在不行可以用fortran做循环，然后单独生成ncl画图所需要的数据。
再不行，可以学一下ncl调用fortran的用法。。那是后话了，我先把手头上的做起来再说~有机会多交流交流