matlab实现对全球不规则投影数据的投影转换

前几个专栏我们讨论了几个不规则的投影转换问题，有需要的可以阅读以下文章：

matlab实现对极地投影数据的投影转换_matlab极地投影-CSDN博客

联合matlab和Arcgis进行netcdf格式的雪覆盖数据的重新投影栅格-CSDN博客

这次遇到的问题是一个墨卡托投影的数据，需要转成标准的WGS经纬度。

使用的数据是HYCOM输出的全球海洋分层温度数据，下载地址：

HYCOM + NCODA Global 1/12° Analysis (GLBa0.08/expt_90.8)

这个数据有一个问题，就是使用的是墨卡托投影，即在南北极变形特别大。见下图：

在读取数据的过程中，我发现其中的经度明显不正常，即有些超过360。

仔细研究后发现，是因为经度需要进行取模处理，即 modulo        = '360 degrees'。

ncdisp显示的变量列表：

Source:
           C:\Users\Hello World!!\DESKTOP\projection\archv.2010_296_00_3zt.nc4
Format:
           netcdf4_classic
Global Attributes:
           _NCProperties = 'version=1|netcdflibversion=4.4.1.1|hdf5libversion=1.8.20'
           Conventions   = 'CF-1.0'
           title         = 'HYCOM GLBa0.08'
           institution   = 'Naval Research Laboratory'
           source        = 'HYCOM archive file'
           experiment    = '90.8'
           history       = 'archv2ncdf3z'
Dimensions:
           MT    = 1     (UNLIMITED)
           Y     = 3298
           X     = 4500
           Depth = 33
Variables:
    MT         
           Size:       1x1
           Dimensions: MT
           Datatype:   double
           Attributes:
                       long_name = 'time'
                       units     = 'days since 1900-12-31 00:00:00'
                       calendar  = 'standard'
                       axis      = 'T'
    Date       
           Size:       1x1
           Dimensions: MT
           Datatype:   double
           Attributes:
                       long_name      = 'date'
                       units          = 'day as %Y%m%d.%f'
                       C_format       = '%13.4f'
                       FORTRAN_format = '(f13.4)'
    Depth      
           Size:       33x1
           Dimensions: Depth
           Datatype:   single
           Attributes:
                       standard_name = 'depth'
                       units         = 'm'
                       positive      = 'down'
                       axis          = 'Z'
    Y          
           Size:       3298x1
           Dimensions: Y
           Datatype:   int32
           Attributes:
                       point_spacing = 'even'
                       axis          = 'Y'
    X          
           Size:       4500x1
           Dimensions: X
           Datatype:   int32
           Attributes:
                       point_spacing = 'even'
                       axis          = 'X'
    Latitude   
           Size:       4500x3298
           Dimensions: X,Y
           Datatype:   single
           Attributes:
                       standard_name = 'latitude'
                       units         = 'degrees_north'
    Longitude  
           Size:       4500x3298
           Dimensions: X,Y
           Datatype:   single
           Attributes:
                       standard_name = 'longitude'
                       units         = 'degrees_east'
                       modulo        = '360 degrees'
    temperature
           Size:       4500x3298x33x1
           Dimensions: X,Y,Depth,MT
           Datatype:   single
           Attributes:
                       coordinates   = 'Longitude Latitude Date'
                       standard_name = 'sea_water_potential_temperature'
                       units         = 'degC'
                       _FillValue    = 1.267650600228229e+30
                       valid_range   = [-7.69922      34.2409]
                       long_name     = '  temp [90.8H]'

因此采用以下的代码进行运算，大于360的，对360度取模，并进行重采样处理，得到规则的经纬度格网。

matlab代码：

file = 'archv.2010_296_00_3zt.nc4';
ncdisp(file)
x = ncread(file,'Longitude');
for i = 1:4500
    for j = 1:3298
        if(x(i,j)>360)
            x(i,j) = mod(x(i,j),360);
        end
    end
end
y = ncread(file,'Latitude');
z = ncread(file,'temperature');
O.lon = x;O.lat = y;O.rg = z(:,:,2);
rg_plot(O)

lonx = reshape(x,4500*3298,1);
laty = reshape(y,4500*3298,1);
rgz  = reshape(O.rg,4500*3298,1);

lon = 0.5:1:359.5;
lat = -89.5:1:89.5;
[lon,lat] = meshgrid(lon,lat);
%% 内插处理
fxy = scatteredInterpolant(double(lonx),double(laty),rgz,'natural');
sla = fxy(lon,lat);
figure
imagesc(sla(:,:,4500))

imagesc(O.rg)

结果图

♥欢迎点赞收藏♥