- 积分
- 86
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-6-27
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 sinohack 于 2015-5-2 18:16 编辑
有CMIP5模式下载的气温nc文件,需要插值求已知站点的温度值,采用反距离加权插值法(IDW)进行插值。经过从论坛和网上学习,用以下下程序进行,根据自己的理解做了注释,但总是在插值的时候出错,请各位帮忙给看看。excel文件为已知站点的经纬度,由两列,一列为经度,一列为纬度。先取了9个站点进行调试程序。
- % 读取nc文件
- ncid = netcdf.open('C:\Users\chris\Desktop\matlab格点数据\tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','NOWRITE'); %打开nc文件
- ncdisp('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc'); %在命令窗中显示nc文件的变量
- tasData = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','tas'); %读入变量tas
- TimeData = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','time'); %读入变量time
- LonData = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','lon'); %读入变量lon
- LatData = ncread('tas_Amon_CCSM4_rcp45_r1i1p1_200601-210012.nc','lat'); %读入变量lat
- [xx,yy]=meshgrid(LonData,LatData);
- tasData=permute(tasData, [2 1 3]);
- %读取nc文件变量完成
- %读取要求站点excl文件经纬度
- lat1 = xlsread('C:\Users\chris\Desktop\matlab格点数据\站点.xls',2,'d2:d10');
- lon1 = xlsread('C:\Users\chris\Desktop\matlab格点数据\站点.xls',2,'e2:e10');
- %读取excl文件经纬度完成
- %IDW(反距离加权插值法)
- %其中x,y,z为已知坐标及其函数值,X,Y为要插值的坐标
- %x,y,z,X,Y最高为二维的,不可为三维
- %不考虑x,y中出现重复坐标的情况
- X=lat1;
- Y=lon1;
- x=xx;
- y=yy;
- z=tasData;
- [m0,n0]=size(x);%返回矩阵x的行数m0,列数n0;[r,c]=size(A)%将矩阵A的行数返回到第一个输出变量r,将矩阵的列数返回到第二个输出变量c
- [m1,n1]=size(X);%返回矩阵X的行数m1,列数n0;
- %生成距离矩阵r(m0*m1*n1,n0)
- for i=1:m1
- for j=1:n1
- 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矩阵的每行中
- end
- end
- %定义插值函数
- for i=1:m1
- for j=1:n1
- if find(r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)==0)
- [m2,n2]=find(r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)==0);
- Z(i,j)=z(m2,n2);
- else
- numerator=sum(sum(z./r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)));
- %z乘距离每一行所有列倒数再竖向求和得到一个行向量,然后再求行向量所有数值的和,r(i,:)表示矩阵r的第i行所有列
- %(sum(x,2)表示矩阵x的横向相加,求每行的和,结果是列向量。而缺省的sum(x)就是竖向相加,求每列的和,结果是行向量。)
- denominator=sum(sum(1./r(m0*n1*(i-1)+m0*(j-1)+1:m0*n1*(i-1)+m0*(j),:)));
- %求距离矩阵倒数,然后竖向相加,得到一个行向量,然后再求行向量各个数的和
- Z(i,j)=numerator/denominator;%IDW(反距离加权插值法)Z=Σ(距离的倒数/距离的倒数之和)*z
- end
- end
- end
复制代码
|
|