ncl的双线性插值问题

afei790912

各位ncl大神们，我有一个ERA-Intrim的相对湿度文件，其信息为：
Variable: r
Type: short
Total Size: 784370400 bytes
            392185200 values
Number of Dimensions: 4
Dimensions and sizes:    [time | 365] x [level | 37] x [latitude | 121] x [longitude | 240]
Coordinates:
            time: [973008..981744]
            level: [1..1000]
            latitude: [90..-90]
            longitude: [ 0..358.5]
Number Of Attributes: 7
  scale_factor :    0.003077844391537395
  add_offset :    79.43179447471586
  _FillValue :    -32767
  missing_value :    -32767
  units :    %
  long_name :    Relative humidity
  standard_name :    relative_humidity
我现在想将之用双线性插值法linint2_points插值到中国区域的探空站点上，并且为每一个站点每一天生成一个txt文件，每一个txt文件为某一站点每一天的37层的相对湿度数据，其ncl代码为：
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"
begin

nf = addfile("~/program/Data/rh.2011.00.nc","r")  ;以只读方式读取文件
lon = nf->longitude(:)    ;将文件中的longtitude坐标变量赋值给变量lon，经度   
latori = nf->latitude(:)  ;将文件中的latitude坐标变量赋值给变量latori，纬度

;将latori变量按递增的方式赋值给变量lat
;将原始纬度数组按从小到大的顺序重新排序
lat = new(121, float)    ;创建一个容量为121，元素类型为float的数组
do i = 0,120
   lat(i) = latori(120-i)
end do

level = nf->level(:)  ;将文件中的levelis变量赋值给变量level
r = short2flt(nf->r(:,:,:,:))       ;将文件中所有范围的r变量赋值给变量r
;读取文本文件中的第1、4、5列赋值给变量station、lon1和lat1
stationPara = "~/program/StationPara.ini"
row = numAsciiRow(stationPara)   ;获取文本文件中的行数

lon1 = new(row,"float")
lat1 = new(row,"float")

obs = asciiread(stationPara,-1,"string")
delim = ","
station = str_get_field(obs,1,delim)   ;用以文件命名
lon1 = stringtofloat(str_get_field(obs,4,delim))
lat1 = stringtofloat(str_get_field(obs,5,delim))

;进行双线性插值计算
newR = linint2_points(lon,lat,r,True,lon1,lat1,0)

;获取nc数据中的时间，转换成年、月、日，用以文件命名
time = nf->time
utc_date = cd_calendar(time,0)
year = tointeger(utc_date(:,0))
month = tointeger(utc_date(:,1))
day = tointeger(utc_date(:,2))

;将newt写成文本文
;按时间、站点输出文本文件
Time = nf->time(:)
TSize = dimsizes(Time)  ;天数
lon1Size =  dimsizes(lon1)  ;站点数
levelSize = dimsizes(level)  ;层数

do i = 0,TSize-1
  do j = 0,lon1Size-1
     interpolationR = newR(i,:,j)
     flipLevel = new(levelSize,"float")   ;存放翻转后level数组   
     flipInterpolationR = new(levelSize,"float")   ;存放翻转后interpolationR数组
     result = new(levelSize,"string")   ;存放数组flipLevel、flipInterpolationR相连的数组
     do k = 0,levelSize-1
        flipLevel(k) = level(levelSize-1-k)
         flipInterpolationR(k) = interpolationR(levelSize-1-k)
        result(k) = sprintf("%8.2f ",flipLevel(k))+sprintf("%8.2f ",flipInterpolationR(k))
     end do

     fName =station(j)+"_"+year(i)+"_"+month(i)+"_"+day(i)+".txt"
     asciiwrite(fName,result)
  end do
end do
end

但插值出来的结果跟探空数据相差特别大，并且连变化趋势都不一致，我怀疑是插值除了问题，但不知如何解决，往各位大神不吝赐教。

并且，还发现一个问题，那就是如果取网格不同的经纬度范围，所取的经纬度范围远远超过所需插值的探空站点经纬度边界，但插值出来的结果不一样，不知怎么回事，往大神指教

longlivehj

r = short2flt(nf->r(:,:,::-1,:))

忽忽儿

非常感谢本帖对于我编格点插值到站点程序的帮助！“r = short2flt(nf->r(:,:,::-1,:))”表示是把r变量所对应的纬度也转换成顺序排列，所以如果不转结果是错误的！虽然能够运行成功数据也不对。感谢大神！

尽头的尽头

请不要重复发帖~~~

尽头的尽头

请不要重复发帖~~~

namelst

值得学习一下

afei790912

尽头的尽头 发表于 2014-9-12 18:58
请不要重复发帖~~~

啊，我重复发帖了吗，那实在是不好意思，

afei790912

longlivehj 发表于 2014-9-12 22:41
r = short2flt(nf->r(:,:,::-1,:))

太感谢了，，我马上试试看，另能否解答一下：所取经纬度范围不同，插值的结果会不一样的问题，如
中国的经纬度范围为：
经度范围：73°33′E至135°05′E
纬度范围：3°51′N至53°33′N
在插值时，ERA-Intrim的经纬度范围我取
经度范围：60°00′E至150°00′E
纬度范围：-15°00′N至70°00′N
该范围已经远远的超过了中国区域，但我取这样的一个格点区域和取全球格点区域对中国的探空站点进行双线性插值，为什么插值出来的结果会不一致呢？
在我的理解中，双线性插值是取站点周边的4个格点值进行插值，那么取这两个范围的格点进行插值的结果应该是一样的，不知怎么回事，麻烦能解答一下，谢谢

longlivehj

afei790912 发表于 2014-9-13 12:56
太感谢了，，我马上试试看，另能否解答一下：所取经纬度范围不同，插值的结果会不一样的问 ...

从你贴出的代码来看，是用全球数据做的，不知道“r = short2flt(nf->r(:,:,::-1,:))”改过之后，结果如何？
至于你说的输入数据范围的问题，结果应该是一样的才对。能贴出你指定输入数据范围时的代码吗？

afei790912

longlivehj 发表于 2014-9-13 16:55
从你贴出的代码来看，是用全球数据做的，不知道“r = short2flt(nf->r(:,:,::-1,:))”改过之后，结果如何 ...

奇怪了，我今天分不同范围又试验了一次，结果居然一样了，，试验情况如下：
估计还是这一句起的作用：“r = short2flt(nf->r(:,:,::-1,:))”

全球取一天的情况：
与上面的代码变化情况：         
time = nf->time(0:1)           

r = short2flt(nf->r(0:1,:,::-1,:))       ;将文件中指定范围的r变量赋值给变量r      

插值结果：中国区域站点50527 2013年1月1日12时的插值结果：
1000.00    91.66
975.00    91.54
950.00    91.67
925.00    90.03
900.00    91.68
875.00    94.38
850.00    94.95
825.00    91.39
800.00    84.23
775.00    82.23
750.00    81.79
700.00    84.10
650.00    87.10
600.00    87.64
550.00    86.18
500.00    88.66
450.00    91.44
400.00    86.36
350.00    96.32
300.00    69.31
250.00    26.84
225.00    17.22
200.00    11.60
175.00     6.91
150.00     5.06
125.00     4.21
100.00     3.91
  70.00     2.52
  50.00     1.15
  30.00     0.07
  20.00     0.01
  10.00     0.01
   7.00     0.00
   5.00     0.00
   3.00     0.00
   2.00     0.00
   1.00     0.00

取包含中国区域一天的情况：
经度范围为：60度~150度
纬度范围为：-15度~70.5度
time = nf->time(0:0)   
lon = nf->longitude(40:100)   ;第40个经度为60度，第100个经度为150度   
                                                                                 
latori = nf->latitude(13:70)  ;第13个纬度为70.5度，第70个纬度为-15度

r = short2flt(nf->r(0:0,:,13:70:-1,40:100))       ;将文件中指定范围的r变量赋值给变量r

插值结果，与上的日期、站点一样：
1000.00    91.66
975.00    91.54
950.00    91.67
925.00    90.03
900.00    91.68
875.00    94.38
850.00    94.95
825.00    91.39
800.00    84.23
775.00    82.23
750.00    81.79
700.00    84.10
650.00    87.10
600.00    87.64
550.00    86.18
500.00    88.66
450.00    91.44
400.00    86.36
350.00    96.32
300.00    69.31
250.00    26.84
225.00    17.22
200.00    11.60
175.00     6.91
150.00     5.06
125.00     4.21
100.00     3.91
  70.00     2.52
  50.00     1.15
  30.00     0.07
  20.00     0.01
  10.00     0.01
   7.00     0.00
   5.00     0.00
   3.00     0.00
   2.00     0.00
   1.00     0.00

afei790912

longlivehj 发表于 2014-9-13 16:55
从你贴出的代码来看，是用全球数据做的，不知道“r = short2flt(nf->r(:,:,::-1,:))”改过之后，结果如何 ...

太感谢您了，

但我感觉这个插值结果仍然是有问题的，比如在450hpa高度附近的相对湿度居然变这么大，没搞懂怎么回事，或许原始数据就有问题吧，如下所示：
1000.00    91.66
975.00    91.54
950.00    91.67
925.00    90.03
900.00    91.68
875.00    94.38
850.00    94.95
825.00    91.39
800.00    84.23
775.00    82.23
750.00    81.79
700.00    84.10
650.00    87.10
600.00    87.64
550.00    86.18
500.00    88.66
450.00    91.44
400.00    86.36
350.00    96.32
300.00    69.31
250.00    26.84
225.00    17.22
200.00    11.60
175.00     6.91
150.00     5.06
125.00     4.21
100.00     3.91
  70.00     2.52
  50.00     1.15
  30.00     0.07
  20.00     0.01
  10.00     0.01
   7.00     0.00
   5.00     0.00
   3.00     0.00
   2.00     0.00
   1.00     0.00