Global mapper 生成的viewshed analysis 的Shapefile格式的文件转换成方位距离格式的文件matlab程序-Go服务器开发

以下是将Global Mapper生成的Viewshed Analysis Shapefile格式文件转换为方位距离格式文件的MATLAB程序：

% 读取Viewshed Analysis Shapefile格式文件
[filename, pathname] = uigetfile('*.shp', 'Select Viewshed Analysis Shapefile');
shpFilePath = fullfile(pathname, filename);
[s, a] = shaperead(shpFilePath);

% 获取DEM文件路径和分辨率信息
demFile = s(1).SOURCE_DEM;
resX = s(1).GRID_DIST_UNITS_METER_X;
resY = s(1).GRID_DIST_UNITS_METER_Y;

% 计算DEM范围和网格大小
xMin = min(a.X(:,1));
xMax = max(a.X(:,1));
yMin = min(a.Y(:,1));
yMax = max(a.Y(:,1));
gridSizeX = round((xMax - xMin) / resX);
gridSizeY = round((yMax - yMin) / resY);

% 创建空的方位距离矩阵
azimuthDistMatrix = zeros(gridSizeY, gridSizeX);

% 遍历每个点并计算其到观测点的方位角和距离
for i=1:length(s)
    obsPointCoord = [s(i).POINT_X, s(i).POINT_Y];
    obsPointElev = s(i).ELEVATION;
    
    % 计算观测点到DEM范围内每个点的距离和方位角
    for j=1:size(azimuthDistMatrix, 2)
        for k=1:size(azimuthDistMatrix, 1)
            currCoord = [xMin + resX * (j - 1), yMin + resY * (k - 1)];
            currElev = demread(demFile, 'coordinate', currCoord);
            if ~isnan(currElev)
                [azimuth, distance] = calcAzimuthDist(obsPointCoord, currCoord);
                azimuthDistMatrix(k,j) = max(azimuthDistMatrix(k,j), distance);
            end
        end
    end
    
    % 显示进度条
    progressText = sprintf('Processing observation point %d of %d...', i, length(s));
    waitbar(i/length(s), hWaitBar, progressText);
end

% 将方位距离矩阵写入文件
[filename, pathname] = uiputfile('*.mat', 'Save As');
if filename ~= 0
    save(fullfile(pathname, filename), 'azimuthDistMatrix');
end

% 计算两点之间的距离和方位角（使用Haversine公式）
function [azimuth, distance] = calcAzimuthDist(coord1, coord2)
    R = 6371e3; % 地球半径，单位：米
    lat1 = deg2rad(coord1(2));
    lat2 = deg2rad(coord2(2));
    deltaLat = deg2rad(coord2(2) - coord1(2));
    deltaLon = deg2rad(coord2(1) - coord1(1));

    a = sin(deltaLat/2)^2 + cos(lat1) * cos(lat2) * sin(deltaLon/2)^2;
    c = 2 * atan2(sqrt(a), sqrt(1-a));

    azimuthRad = atan(sin(deltaLon) / (cos(lat1) * tan(lat2) - sin(lat1) * cos(deltaLon)));
    if (lat2 < lat1)
        azimuthRad = pi + azimuthRad;
    end

    azimuth = rad2deg(azimuthRad);
    distance = R * c;
end

使用方法：