matlab实现对极地投影数据的投影转换

news2024/11/17 0:21:09

图片摘要

前段时间，我写了一个专栏《联合matlab和Arcgis进行netcdf格式的雪覆盖数据的重新投影栅格》，介绍了如何利用Arcgis实现极地投影的转换。今天，我提供一个matlab重采样的方法实现投影转换。

我们这次使用的数据是北冰洋逐日的海平面高度数据。但是我们发现其是极地投影，非常不利用后续的计算分析。我也在arcgis中尝试过，但是数据存在的区域远远大于实际合理的经纬度范围。如下图。真实际的范围应该是与中间的圆形边界相切。

于是，在matlab中，我们可以首先读取数据，然后得到内部相切的圆形区域。

file = ('dt_arctic_multimission_v1.1_sea_level_20160701_20190429.nc');
ncdisp(file)
lon = ncread(file,'longitude');
lat = ncread(file,'latitude');
sla1 = ncread(file,'sla');
time = ncread(file,'time');
P.lon = lon;P.lat = lat;
P.rg  = sla1(:,:,1);
imagesc(P.rg)
%% time convert
dt = datetime((time)*24*3600, 'ConvertFrom', 'epochtime', 'Epoch', '1950-01-01');
[y,m,d] = ymd(dt);

xv = -179.8358:0.25:179.8358;
yv = 50.0012:0.25:89.8417;
[xv1,yv1] = meshgrid(xv,yv);

lon1 = lon(186:535,186:535);
lat1 = lat(186:535,186:535);
lat1(lat1>1000)=0;
lon1(lon1>1000)=0;
xx = reshape(lon1,350*350,1);
yy = reshape(lat1,350*350,1);

接下来是采用内插的方法，得到规则化格网对应的数据。

N = sla1(:,:,1);
N(N>1000)=0;
%% 得到与圆形边界相切的区域
N1  = N(186:535,186:535);
zz = reshape(N1,350*350,1);
out = [xx,yy,zz];
ind = find(out(:,1)==0);
%% 将没有意义的经纬度格网赋值为0
out(ind,:)=[];
out = double(out);
%% 内插处理
fxy = scatteredInterpolant(out(:,1),out(:,2),out(:,3),'natural');
sla(:,:,i) = fxy(xv1,yv1);