MATLAB读取nc文件，并绘制南海水深图

记录第一次拿MATLAB画出能看的图！

---

1. 代码实现

clear all
clc
%%读取数据
Inpath = 'D:\MATLAB\bathymetry\ETOPO_2022_v1_60s_N90W180_bed.nc'
ncdisp(Inpath)
lon = ncread(Inpath,'lon');
lat = ncread(Inpath,'lat');
DEPTH = ncread(Inpath,'z');
lon_real = lon(find(lon>= 109 & lon<= 123));
lat_real = lat(find(lat >= 11 & lat <= 23));
DEPTH(find(DEPTH > 0))=nan;
DEPTH_real=DEPTH(find(lon>= 109 & lon<= 123),find(lat >= 11 & lat <= 23));
%%可视化
DEPTH_real=DEPTH_real';
[LON,LAT] = meshgrid(lon_real,lat_real);
hold on;
figure(1)
set(0,'defaultfigurecolor','w')   %%% 图片底色为白色
m_proj('lambert','lon',[109  123],'lat',[11 23]);  %%% Miller 投影
hold on
%%m_plot(108.58337402,18.18760655,'.','MarkerSize',20,'MarkerEdgeColor','r')
m_contourf(lon_real,lat_real,DEPTH_real)
m_gshhs_f('patch',[.7 .7 .7],'edgecolor','none');
m_coast('patch',[.6 .6 .6],'edgecolor','k'); 
m_grid('linest','none','tickdir','out','box','fancy','fontsize',16) 
colorbar
colormap(m_colmap('blues'));  
set(gcf,'color','w');

---

2. 读取数据

ncdisp(‘filepath’)：显示nc文件中变量维度、变量名称

lon = ncread(Inpath,'lon');
lat = ncread(Inpath,'lat');
DEPTH = ncread(Inpath,'z');

从nc文件中读取经纬度、深度数据，并以矩阵的形式赋值给变量lon,lat,DEPTH。

注意：这里的赋值语句后面要加上分号，作用是不显示赋值结果；否则经纬度数据太多，下面的命令行窗口要运行很久。

lon_real = lon(find(lon>= 109 & lon<= 123));
lat_real = lat(find(lat >= 11 & lat <= 23));

所下载的nc文件中的经纬度数据是全球范围的，我只需要南海的经纬度数据，因此用find（）函数提取满足条件的元素下标。例如，find(lon>= 109 & lon<= 123)，运行结果是：找到满足该条件的经度下标，并把下标储存在ans中。用lon(find(lon>= 109 & lon<= 123))将ans中储存的下标还原成元素值，就得到了我需要的经度：lon_real。

 DEPTH(find(DEPTH > 0))=nan；

高程大于0的数据是陆地数据，属于冗余数据，把他们的值变成nan，这样程序运行时会跳过这些数据。

DEPTH_real=DEPTH(find(lon>= 109 & lon<= 123),find(lat >= 11 & lat <= 23))；

水深数据是二维矩阵，它的行和列刚好和经度、维度数据的下标对应，因此通过这行代码可以提取出南海范围的水深数据。

 [LON,LAT] = meshgrid(lon_real,lat_real);

3. 画水深图

figure（1）：打开一个画图窗口，程序显示结果如下图：

 set(0,'defaultfigurecolor','w') ：设置图片的背景颜色为白色

m_proj('lambert','lon',[109  123],'lat',[11 23]);

设置投影方式为Lambert投影，设置投影的经纬度范围。

m_contourf(lon_real,lat_real,DEPTH_real)     %%最重要的一步！
contourf(X,Y,Z)函数：绘制矩阵Z的等值线图，并且指定矩阵Z中各值的x,y坐标。

 m_coast('patch',[.6 .6 .6],'edgecolor','k'); 

设置海岸线，‘patch’代表填充，[.6 .6 .6]控制填充颜色，三个RGB参数分别是红绿蓝。

'edgecolor','k':设置岸线轮廓线为黑色。这里可以看到岸线轮廓非常粗糙，因此正在学习使用m_gshhs_f（）函数来绘制精细岸线。下一篇文章中已解决。

 m_grid('linest','none','tickdir','out','box','fancy','fontsize',16)

生成经纬度网格，对其中参数的说明:

linest：网格线型，有-- -. : - 四种

tickdir：刻度线方向。对于fancy来说in表示黑白框out表示线条框

box是边框设置：on有，off无，fancy为黑白框

colorbar
显示色阶的颜色栏

colormap(m_colmap('blues'));  
将颜色栏设置为蓝色，按颜色深浅区分

set(gcf,'color','w');

 例如set(gcf,'color','k')执行结果：

---

附录：数据下载教程

水深文件下载网址：ETOPO Global Relief Model | National Centers for Environmental Information (NCEI) (noaa.gov)https://www.ncei.noaa.gov/products/etopo-global-relief-model#:~:text=The%20ETOPO%20Global%20Relief%20Model%20integrates%20topography%2C%20bathymetry%2C,renderings%20of%20geophysical%20characteristics%20of%20the%20earth%E2%80%99s%20surface.

进入该网站后向下拉，找到这一部分：

15/30/60 Arc-Second Resolution代表不同的分辨率，Arc-Second即经纬度单位中，度分秒中的“秒”，因此15 Arc-Second Resolution是最精细的数据。本文中使用的数据为：Bedrock-elevation netCDF。