关于nc数据插值到站点的问题

sinohack

有CMIP5模式下载的气温nc文件，需要插值求已知站点的温度值，采用反距离加权插值法（IDW）进行插值。经过从论坛和网上学习，用以下下程序进行，根据自己的理解做了注释，但总是在插值的时候出错，请各位帮忙给看看。excel文件为已知站点的经纬度，由两列，一列为经度，一列为纬度。先取了9个站点进行调试程序。

1. 
   
2. % 读取nc文件
   
3. ncid = netcdf.open('C:\Users\chris\Desktop\matlab格点数据\tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','NOWRITE'); %打开nc文件
   
4. ncdisp('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc'); %在命令窗中显示nc文件的变量
   
5. tasData  = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','tas'); %读入变量tas
   
6. TimeData  = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','time'); %读入变量time
   
7. LonData  = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','lon'); %读入变量lon
   
8. LatData  = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','lat'); %读入变量lat
   
9. [xx,yy]=meshgrid(LonData,LatData);
   
10. tasData=permute(tasData, [2 1 3]);
    
11. %读取nc文件变量完成
    
12. 
    
13. 
    
14. %读取要求站点excl文件经纬度
    
15. lat1 = xlsread('C:\Users\chris\Desktop\matlab格点数据\站点.xls',2,'d2:d10');
    
16. lon1 = xlsread('C:\Users\chris\Desktop\matlab格点数据\站点.xls',2,'e2:e10');
    
17. %读取excl文件经纬度完成
    
18. 
    
19. 
    
20. %IDW（反距离加权插值法）
    
21. %其中x,y,z为已知坐标及其函数值,X,Y为要插值的坐标
    
22. %x,y,z,X,Y最高为二维的，不可为三维
    
23. %不考虑x，y中出现重复坐标的情况
    
24. X=lat1;
    
25. Y=lon1;
    
26. x=xx;
    
27. y=yy;
    
28. z=tasData;
    
29. 
    
30. 
    
31. 
    
32. [m0,n0]=size(x);%返回矩阵x的行数m0，列数n0；[r,c]=size（A）%将矩阵A的行数返回到第一个输出变量r，将矩阵的列数返回到第二个输出变量c
    
33. [m1,n1]=size(X);%返回矩阵X的行数m1，列数n0；
    
34. %生成距离矩阵r(m0*m1*n1,n0)
    
35. for i=1:m1
    
36.     for j=1:n1
    
37.         r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)=sqrt((X(i,j)-x).^2+(Y(i,j)-y).^2);%求各点的距离并存入r矩阵的每行中
    
38.     end
    
39. end
    
40. %定义插值函数
    
41. for i=1:m1
    
42.     for j=1:n1
    
43.         if find(r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)==0)
    
44.             [m2,n2]=find(r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)==0);
    
45.             Z(i,j)=z(m2,n2);
    
46.         else
    
47.             numerator=sum(sum(z./r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)));
    
48.             %z乘距离每一行所有列倒数再竖向求和得到一个行向量，然后再求行向量所有数值的和，r(i,:)表示矩阵r的第i行所有列
    
49.             %（sum(x,2)表示矩阵x的横向相加，求每行的和，结果是列向量。而缺省的sum(x)就是竖向相加，求每列的和，结果是行向量。）
    
50.             denominator=sum(sum(1./r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)));
    
51.             %求距离矩阵倒数，然后竖向相加，得到一个行向量，然后再求行向量各个数的和
    
52.             Z(i,j)=numerator/denominator;%IDW（反距离加权插值法）Z=Σ（距离的倒数/距离的倒数之和）*z
    
53.         end
    
54.     end
    
55. end
    

复制代码

歌留多

楼主，我想问一下，你最后插值得到的Z的矩阵的规格是怎样的，跟X一样么？

鱼香Jack

~~~~~~~~~~~~~~~~~~~~~~~~~~

coldearth

学习了，非常感谢！

kimhyunjoong

请问楼主最后是如何解决的？
我运行时也是出现错误提示：

1. 错误使用  ./
   
2. 矩阵维度必须一致。
   
3. 
   
4. 出错 chazhi (line 47)
   
5.             numerator=sum(sum(z./r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)));

复制代码

ztf6863

谢谢楼主分享啦

smilelxf

请问一下，如何把CMIP模式插值到2.5°*2.5*的分辨率上

loser.

kimhyunjoong 发表于 2017-8-27 16:06
请问楼主最后是如何解决的？
我运行时也是出现错误提示：

请问，您解决这个问题了么

kimhyunjoong

loser. 发表于 2019-6-30 23:30
请问，您解决这个问题了么

时间太久，我也不记得了

loser.

kimhyunjoong 发表于 2019-7-8 14:35
时间太久，我也不记得了

{:5_275:}{:5_275:}