java使用gdal可视域分析，判断返回该坐标是否可视

gdal自带的ViewshedMode默认为一圈圆形通视分析，感觉应该不能实现自定义扇形范围通视分析，因此这里采用射线法实现了简单的扇形通视分析

这控制采样点越多，旋转角越小则分析结果越多越准

package com.leador.gsp.module.parse.test;

import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.gdal;

/**
 * 考虑高度及左右范围可视范围
 */
public class RasterGdalTest_range {
    public static void main(String[] args) {
        double startLon = 4.189;
        double startLat = 22.061;
        double startHeight = 100;
        double directionStart = 30.0;//其实角度
        double directionEnd = 60.0;//终止角度
        double testAngle = directionStart;
        double testDirectionAngle = Math.toRadians(testAngle);
        double distanceMeters = 300000; // 要移动的距离（单位：米）
        double earthRadius = 6371.393 * 1000; // 地球的平均半径（单位：米）
        double minHeight = 0;
        // 计算最小可见高度
        if ((startHeight * startHeight - distanceMeters * distanceMeters) > 0){
            minHeight = startHeight - Math.sqrt(startHeight * startHeight - distanceMeters * distanceMeters);
        }
        // 计算最大可见高度
        double maxHeight = startHeight + Math.sqrt(startHeight * startHeight + distanceMeters * distanceMeters);
        System.out.println("最小可见高度：" + minHeight + " 米");
        System.out.println("最大可见高度：" + maxHeight + " 米");
        while(testAngle >= directionStart && testAngle <= directionEnd){
            testDirectionAngle = Math.toRadians(testAngle);
            // 计算目标点的经度和纬度
            double targetLat = startLat + distanceMeters*Math.cos(testDirectionAngle)/(earthRadius*2*Math.PI/360);
            double targetLon = startLon + distanceMeters * Math.sin(testDirectionAngle)/(earthRadius*Math.cos(targetLat *Math.PI/180 )*2*Math.PI/360);
            System.out.println("目标点经度:  " + targetLon + ", 目标点纬度:  " + targetLat);
            visibleResult(startLon, startLat, targetLon, targetLat, minHeight, maxHeight);
            System.out.println("-----------------");
            //每次角度转五度
            testAngle  =  testAngle + 5;
            //可视高度归零
            minHeight = 0;
        }
    }
    public static void visibleResult (double startLon, double startLat, double targetLon, double targetLat, double minHeight, double maxHeight){
        gdal.AllRegister();  // 注册所有的GDAL驱动
        // 打开DEM数据集
        Dataset dataset = gdal.Open("C:\\Users\\82305\\Desktop\\1\\gebco.tif");
        // 计算直线上的采样点（简化为每隔一定距离取一个采样点）
        int numPoints = 20;
        double[] xCoords = new double[numPoints];
        double[] yCoords = new double[numPoints];
        for (int i = 0; i < numPoints; i++) {
            xCoords[i] = startLon + (targetLon - startLon) / (numPoints - 1) * i;
            yCoords[i] = startLat + (targetLat - startLat) / (numPoints - 1) * i;
        }
        // 获取DEM数据集的地理转换信息
        double[] geotransform = dataset.GetGeoTransform();
        double pixelWidth = geotransform[1];
        double pixelHeight = geotransform[5];
        // 获取DEM数据集的高程信息
        Band band = dataset.GetRasterBand(1);
        int[] elevations = new int[numPoints];
        for (int i = 0; i < numPoints; i++) {
            //TODO （采样点的x - 地理左上角坐标x）/宽度  得到列索引
            int x = (int) ((xCoords[i] - geotransform[0]) / pixelWidth);
            int y = (int) ((yCoords[i] - geotransform[3]) / pixelHeight);
            int[] sample = new int[1];
            band.ReadRaster(x, y, 1, 1, sample);
            elevations[i] = sample[0];
            //如果碰到第一个采样点高度大于观测点高度即可视 ，后续大于上一个采样点高度也为可视
            if (elevations[i] >= minHeight && elevations[i] < maxHeight){
                System.out.println("可视Point " + i + ": (" + xCoords[i] + ", " + yCoords[i] + ") - Elevation: " + elevations[i]);
                //当看到新的采样点高度后覆盖原有的最小可视高度
                minHeight = elevations[i];
            }else {
                //采样点高度低于观测点高度即不可视
                System.out.println("不可视Point " + i + ": (" + xCoords[i] + ", " + yCoords[i] + ") - Elevation: " + elevations[i]);
            }
        }
        // 关闭DEM数据集
        dataset.delete();
    }
}