m_map导入本地地形数据

m_map绘制地形图时，虽然自带有1°的地形图以及从NOAA下载的1分的地形图（详见：Matlab下地形图绘图包m_map安装与使用），但有时需要对地形图分辨率的要求更高，便无法满足。

此时，需要导入本地地形数据文件，例如多波束测深。

这里以从https://www.ncei.noaa.gov/maps/grid-extract/下载的15秒地形数据为例，展示如何导入本地地形数据并绘图。

关键函数为m_shadedrelief。

1. 1分的分辨率地图

m_proj('mercator','long',[118 124],'lat',[20 24]);

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

m_etopo2('shadedrelief','lightangle',45,'gradient',5);

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

set(gcf,'position',[10 10 1200 800])

figname='test1'
print('-dpng','-r1000',[figname,'.png'])  % 导出png图片

结果如下：

2. 15秒的分辨率地图

使用geotiffread函数读取tif格式的高程数据与坐标。
还需要m_image函数用以校正投影坐标，和m_shadedrelief配合使用。

clear;close all;clc

file='image.tif';
[elev,information]=geotiffread(file);

lonlim=information.LongitudeLimits;
latlim=information.LatitudeLimits;

elev=double(flipud(elev));
[nlat,nlon]=size(elev);
Lon=[linspace(lonlim(1),lonlim(2),nlon)];
Lat=[linspace(latlim(1),latlim(2),nlat)]';

m_proj('mercator','long',[lonlim(1) lonlim(2)],...
    'lat',[latlim(1) latlim(2)]);
caxis([-6000 6000]) 
colormap(slanCM('terrain')) 
hc=colorbar;
set(get(hc,'title'),'string','Elevation(m)')	
set(hc,'tickdir','out','linewid',2)
[IM,X,Y]=m_image(Lon,Lat,elev); % 此命令必须放在m_shadedrelief之前用以校正坐标
m_shadedrelief(X,Y,IM,'coords','map','lightangle',45,'gradient',5)

m_gshhs('ic','color','k')

m_grid('box','on','tickdir','out','gridlines','no','linewid',2)

set(gcf,'position',[10 10 1200 800])

figname='test2'
print('-dpng','-r600',[figname,'.png'])  % 导出png图片

• 可见15秒分辨率地图比1分地图能显示更多细节。

注：绘图色标参考自：https://blog.csdn.net/slandarer/article/details/127719784