【MATLAB】自学记录之读取DEM高程数据文件并渲染成三维地形图

1. 前言

近日在学习MATLAB编程以及地理高程数据处理等相关知识时，希望通过MATLAB的绘图等相关函数，读取高程数据文件，最后以可视化的方式展示全球陆地范围内的三维高程数据图。

2. 运行环境及数据

序号	配置项	说明
1	CPU	Intel i5-12490F
2	内存	16G*2, 3600MHz
3	磁盘	256G，SSD
4	显卡	RTX4060Ti
5	MATLAB	2023b
6	高程数据文件	GMTED2010.TIFF，包含全球陆地的高程数据

3. 实现过程

1. 通过readgeoraster函数读取tiff文件，获取DEM数据大小。

   clc;clear;
   file = 'GMTED2010.tif';
   [dem_data, ~] = readgeoraster(file);
   [row,col,~] = size(dem_data);
   

2. （可选）设置缩放倍数。

   由于计算机内存受限，以原始DEM数据进行渲染3D效果时，会导致内存不够。

   % 设置缩放倍数
   % 10代表从dem_data数据中每隔10个点取一次
   multiple = 10;
   % 缩放后的XY坐标向量（经度轴和纬度轴）
   small_col = col / multiple;
   small_row = row / multiple;
   

3. 初始化Z轴向量（高程轴）。

   % 方法一
   alt = int16(zeros(small_row, small_col));
   
   for ri = 1: small_row
       for ci= 1: small_col
           alt(ri, ci) = dem_data((ri-1)*multiple + 1, (ci-1)*multiple + 1);
       end
       disp("total data: " + small_row*small_col + ", copied: " + ri*ci);
   end
   
   % 方法二
   alt = dem_data(1:multiple:end,1:multiple:end, 1);
   

4. （可选）清除无用的变量，腾出内存空间。

   clear dem_data;
   

5. 对横坐标longitude和纵坐标latitude生成二维网格矩阵，矩阵大小为 $small\_row \times small\_col$，与alt变量维度保持一致。

   [longitude, latitude] = meshgrid(1:small_col, 1:small_row);
   
   % 默认是double类型，转换成int16，减少内存占用
   longitude = int16(longitude);
   latitude = int16(latitude);
   

6. 使用mesh函数渲染3D曲面图形。

   % 将alt矩阵进行左右翻转
   alt = fliplr(alt);
   mesh(longitude, latitude, alt);
   

7. 为图形加上三轴标签。

   xlabel("longitude");
   ylabel("latitude");
   zlabel("alt");
   hold on;
   

4. 渲染效果

5. 完整代码

clc;clear;
file = 'GMTED2010.tif';
[dem_data, ~] = readgeoraster(file);
[row,col,~] = size(dem_data);

multiple = 10;

small_col = col / multiple;
small_row = row / multiple;

% alt = int16(zeros(small_row, small_col));

%for ri = 1: small_row
%    for ci= 1: small_col
%        alt(ri, ci) = dem_data((ri-1)*multiple + 1, (ci-1)*multiple + 1);
%    end
%    disp("total data: " + small_row*small_col + ", copied: " + ri*ci);
%end

alt = dem_data(1:multiple:end, 1:multiple:end, 1);

clear dem_data;
[longitude, latitude] = meshgrid(1:small_col, 1:small_row);

longitude = int16(longitude);
latitude = int16(latitude);
alt = fliplr(alt);
mesh(longitude, latitude, alt);
xlabel("longitude");
ylabel("latitude");
zlabel("alt");
hold on;