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)