m_map绘制多波束数据

之前，分享了m_map如何绘制本地高程数据：m_map导入本地地形数据。

其实，m_map还可以绘制多波束测深数据。

1. 多波束网站

多波束测深数据可以从https://www.ncei.noaa.gov/maps/grid-extract/或者https://www.ncei.noaa.gov/maps/bathymetry/下载。
多波束测深数据的分辨率为3秒，约为90 m 。

在这里插入图片描述

2. 绘图

以纽约以东的一片海域为例。

这是原始的多波束数据，里面无有效数据处，被填充为0值，显示如下：

在这里插入图片描述

将0值替换为NaN值，即显示为空白，显示如下：

在这里插入图片描述

将m_map自带的1分的分辨率地形图，融合到多波束数据中，显示如下：

在这里插入图片描述

显示带有海山的区域：

4. 代码

mat_name='newyork';
fignum=1

load([mat_name,'_elevation.mat'])
elev=flipud(elev);

% 将0值替换为NaN
% elev(find(elev==0))=nan;

lonlim=[extent(1:2)]
latlim=[extent(3:4)]
[nlat,nlon]=size(elev);
Lon=[linspace(lonlim(1),lonlim(2),nlon)];
Lat=[linspace(latlim(1),latlim(2),nlat)]';


%% 高程数据融合
% 在多波束数据无数据（标志为0处，替换为自带的1分高程数据）
% 首先提取1分的高程数据

[ELEV,LON,LAT]=m_etopo2([lonlim(1)-1 lonlim(2)+1 latlim(1)-1,latlim(2)+1]);
[Lon,Lat]=meshgrid(Lon,Lat);
% 将1分数据插值为和多波束数据同样的大小
ELEV_p=interp2(LON,LAT,ELEV,Lon,Lat);

for i=1:nlat
    for j=1:nlon
        if elev(i,j)==0
            elev(i,j)=ELEV_p(i,j);
        end
    end
end

%% 绘图

m_proj('mercator','long',[-67 lonlim(2) ],'lat',[38.5,latlim(2)]);

caxis([-5000 0])
colormap(slanCM('rainbow')) 
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')
set(hc,'tickdir','out')

m_shadedrelief(Lon(1,:),Lat(:,1),elev,'lightangle',45,'gradient',10) 

m_gshhs('ic','color','k')
m_grid('box','on','tickdir','out','gridlines','no')

set(gcf,'position',[10 10 400 400])

figname=[mat_name,num2str(fignum)];
print('-dpng','-r500',[figname,'.png'])  % 导出png图片