请求的 43200x64800 (20.9GB)数组超过预设的最大数组大小

我的目标提取火山轮廓,其中一步是根据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
  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 7
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值