我的目标提取火山轮廓,其中一步是根据DEM,将火山轮廓理想化成一个圆锥,使圆锥有DEM的坐标系统,并且生成圆锥的Aspect图,半径和高程是根据DEM获取的。裁剪的影像小的时候能运行出来,但是影像大了一点就报错,请问具体要怎么改?
出错提示:
请求的 43200x64800 (20.9GB)数组超过预设的最大数组大小。创建大于此限制的数组可能需要较长时间,并且会导致 MATLAB 无响应。
出错 Untitled2_2>createConeWithDEM (第 15 行)
dem = double(dem);
出错 Untitled2_2 (第 10 行)
createConeWithDEM(dem_file, r, h, n, num_segments, new_dem_file, new_aspect_file);
DEM在这:链接:https://pan.baidu.com/s/1RmzktqFV6HvtZIxIcqFYZg
提取码:1111
clc;
clear all;
dem_file = 'D:\数据\1\extract_dtm.tif';
r = 2300;
h = 331;
n = 50;
num_segments = 10;
new_dem_file = 'D:\数据\1\Cone_DEM1.tif';
new_aspect_file = 'D:\数据\1\Cone_Aspect1.tif';
createConeWithDEM(dem_file, r, h, n, num_segments, new_dem_file, new_aspect_file);
function createConeWithDEM(dem_file, r, h, n, num_segments, new_dem_file, new_aspect_file)
% 读取DEM数据
[dem, R] = geotiffread(dem_file);
dem = double(dem);
[rows, cols] = size(dem);
% 获取原始DEM文件的GeoTIFFTags
dem_info = geotiffinfo(dem_file);
geo_key_directory_tag = dem_info.GeoTIFFTags.GeoKeyDirectoryTag;
% 创建新的DEM文件和坡向文件
cone_dem = zeros(rows, cols);
aspect_dem = zeros(rows, cols);
% 计算圆锥体顶点坐标
cone_x = round(cols / 2);
cone_y = round(rows / 2);
% 循环遍历DEM的每个像素点
for row = 1:rows
for col = 1:cols
% 计算当前像素点到圆锥体顶点的距离
dist = sqrt((row - cone_y)^2 + (col - cone_x)^2);
% 如果当前像素点在圆锥体范围内
if dist <= r
% 计算当前像素点的圆锥体高程
cone_height = h * (1 - dist / r);
% 将圆锥体高程添加到新的DEM文件中
cone_dem(row, col) = cone_height;
% 计算坡向并添加到坡向文件中
dx = col - cone_x;
dy = row - cone_y;
aspect = atan2(dy, dx) * 180 / pi;
aspect = mod(-aspect + 90, 360); % 转换为与正北方向一致的坡向
aspect_dem(row, col) = aspect;
else
% 将原始DEM高程添加到新的DEM文件中
cone_dem(row, col) = dem(row, col);
aspect_dem(row, col) = NaN;
end
end
end
% 保存新的DEM文件和坡向文件,同时保留地理参考信息
cone_dem = uint16(cone_dem);
geotiffwrite(new_dem_file, cone_dem, R, 'GeoKeyDirectoryTag', geo_key_directory_tag);
geotiffwrite(new_aspect_file, aspect_dem, R, 'GeoKeyDirectoryTag', geo_key_directory_tag);
% 绘制3D圆锥体(原代码部分)
theta = linspace(0, 2 * pi, n);
z_segments = linspace(0, h, num_segments);
vertices = [];
colors = [];
for z = z_segments
[x_temp, y_temp] = pol2cart(theta, r * (1 - z / h));
% 考虑DEM的空间位置
x_temp = x_temp + cone_x;
y_temp = y_temp + cone_y;
vertices = [vertices; x_temp' y_temp' repmat(z, length(x_temp), 1)];
height_normalized = z / h;
cmap = parula;
color = cmap(round(height_normalized * (size(cmap, 1) - 1)) + 1, :);
colors = [colors; repmat(color, length(x_temp), 1)];
end
% 添加顶点
vertices = [vertices; cone_x cone_y h];
colors = [colors; cmap(end, :)];
% 计算正圆锥体的面片索引
tri = [];
n_segments = length(z_segments);
for i = 1:n_segments - 1
for j = 1:n - 1
tri = [tri; n * (i - 1) + j, n * (i - 1) + j + 1, n * i + j];
tri = [tri; n * i + j, n * (i - 1) + j + 1, n * i + j + 1];
end
tri = [tri; n * (i - 1) + n, n * (i - 1) + 1, n * i + n];
tri = [tri; n * i + n, n * (i - 1) + 1, n * i + 1];
end
for j = 1:n - 1
tri = [tri; n * (n_segments - 1) + j, n * (n_segments - 1) + j + 1, n * n_segments + 1];
end
tri = [tri; n * (n_segments - 1) + n, n * (n_segments - 1) + 1, n * n_segments + 1];
% 绘制3D正圆锥体
trisurf(tri, vertices(:, 1), vertices(:, 2), vertices(:, 3), 'FaceColor', 'interp', 'FaceVertexCData', colors, 'EdgeColor', 'none');
axis equal;
xlabel('X轴');
ylabel('Y轴');
zlabel('Z轴');
title('3D正圆锥体按高度分层着色');
colormap(parula);
colorbar;
end