GDAL使用

使用gdal.netcore来读取和生成栅格文件。
优点：自带gdal运行时相关文件，不用额外再安装gdal库
缺点：导致发布的文件变大很多，比如Win+Linux的运行时加起来就超过了400M，所以最好是按需加载对应的运行时。比如DEBUG是运行在win的，就添加MaxRev.Gdal.WindowsRuntime.Minimal包，然后Release环境是运行在linux x64的，就添加MaxRev.Gdal.LinuxRuntime.Minimal.x64。

如图所示，仅支持64位系统。
还有就是太老的系统不支持，比如centOS7,那就只能是源码编译安装Gdal(比较折腾的过程)，然后再单独修改gdal注册相关的代码。不过老系统更推荐的做法是使用docker，基础容器镜像使用较新的系统版本，比如debian 12之类的。

可以通过msbuild condition比较编译时的提供的环境变量来引用不同的运行时包

  <ItemGroup Condition="'$(TargetRuntime)'=='win'">
    <PackageReference Include="MaxRev.Gdal.WindowsRuntime.Minimal" />
  </ItemGroup>
  <ItemGroup Condition="'$(TargetRuntime)'=='lin-x64'">
    <PackageReference Include="MaxRev.Gdal.LinuxRuntime.Minimal.x64" />
  </ItemGroup>

然后在dotnet build或publish时传入对应的参数 -p:TargetRuntime=win 或-p:TargetRuntime=lin-x64

计算FVC

计算方法如下，读取卫星遥感影像的3、4波段，计算NDVI值中的5%和95%值作为极值，然后再根据极值计算FVC。
由于单个tif文件可能极大，所以第一次计算NDVI极值时不存计算出的NDVI值（一个省的10m分辨率影像就近100G了）
后续计算FVC时每次只读取缓冲区大小的数据进行计算。
缺点是串行计算十分缓慢，直接使用Task.WhenAll又会导致读取栅格报错protected memory。可能在读取和写入栅格数据时进行加锁处理是可行的。

参数是遥感影像栅格和FVC结果栅格的路径，返回值是遥感影像栅格的大小、坐标系、地理变换等信息。

    /// <summary>
    /// 基于哨兵四波段遥感卫星影像（B\G\R\NIR）,计算植被归一化指数（NDVI） 基于NDVI，计算植被覆盖度FVC
    /// </summary>
    /// <param name="tifPath">遥感影像路径</param>
    /// <param name="resultTifPath">FVC栅格路径</param>
    private (double[] transform, int xSize, int ySize, SpatialReference srs, string projection) Fvc(string tifPath,
        string resultTifPath)
    {
        Console.WriteLine($"{DateTime.Now: mm:ss} 开始计算FVC");
        using var ds = Gdal.Open(tifPath, Access.GA_ReadOnly);
        using var driver = Gdal.GetDriverByName("GTiff");
        using var bandR = ds.GetRasterBand(3);
        using var bandNIR = ds.GetRasterBand(4);
        //NDVI=float(band4-band3)/(band4+band3)
        //FVC=（NDVI-bandmin）/(bandmax-bandmin) bandmin取NDVI中累计值为最接近5%的值，bandmax取NDVI中累计值为最接近95%的值,FVC计算结果应为0~1区间，计算结果中存在小于0的值归为0，存在大于1的值归为1
        //设置缓冲区2048*2048，每次只从4个波段中读取对应缓冲区的数据进行计算
        var xSize = ds.RasterXSize;
        var ySize = ds.RasterYSize;
        var blockSize = 2048;
        var xCount = (int)Math.Ceiling((float)xSize / blockSize);
        var yCount = (int)Math.Ceiling((float)ySize / blockSize);
        bandR.GetNoDataValue(out var noDataValue, out _);
        if (File.Exists(resultTifPath)) File.Delete(resultTifPath);
        var transform = new double[6];
        ds.GetGeoTransform(transform);
        ConcurrentDictionary<double, int> nvdiSet = new();
        for (int y = 0; y < yCount; y++)
        {
            var bufferR = new List<double[]>(xCount);
            var bufferNIR = new List<double[]>(xCount);
            for (int x = 0; x < xCount; x++)
            {
                var offsetX = x * blockSize;
                var offsetY = y * blockSize;
                var actualWidth = Math.Min(blockSize, Math.Max(xSize - offsetX, 0));
                var actualHeight = Math.Min(blockSize, Math.Max(ySize - offsetY, 0));
                bufferR.Add(new double[actualWidth * actualHeight]);
                bandR.ReadRaster(offsetX, offsetY, actualWidth, actualHeight, bufferR[x], actualWidth, actualHeight, 0,
                    0);
                bufferNIR.Add(new double[actualWidth * actualHeight]);
                bandNIR.ReadRaster(offsetX, offsetY, actualWidth, actualHeight, bufferNIR[x], actualWidth, actualHeight,
                    0, 0);
            }

            //计算NDVI
            for (var i = 0; i < xCount; i++)
            {
                var index = i;
                for (var j = 0; j < bufferR[index].Length; j++)
                {
                    var r = bufferR[index][j];
                    var nir = bufferNIR[index][j];
                    if (double.IsNaN(r) || double.IsNaN(nir)) continue;
                    var ndvi = (nir - r) / (nir + r);
                    if (double.IsNaN(ndvi)) continue;
                    nvdiSet.AddOrUpdate(Math.Round(ndvi, 5), 1, (key, oldValue) => oldValue + 1);
                }
            }
        }

        //计算5%和95%对应的Count
        var totalCount = nvdiSet.Sum(d => d.Value);
        var ndviMin = GetPercentCount(totalCount, 0.05, nvdiSet.OrderBy(d => d.Key));
        var ndviMax = GetPercentCount(totalCount, 0.05, nvdiSet.OrderByDescending(d => d.Key));

        using var dsResult = driver.Create(resultTifPath, xSize, ySize, 1, DataType.GDT_Float32, null);
        var projection = ds.GetProjection();
        dsResult.SetProjection(projection);
        dsResult.SetGeoTransform(transform);
        var srs = ds.GetSpatialRef();
        dsResult.SetSpatialRef(srs);
        using var bandResult = dsResult.GetRasterBand(1);
        bandResult.SetNoDataValue(noDataValue);

        var ndviDiff = ndviMax - ndviMin;
        for (int y = 0; y < yCount; y++)
        {
            //计算NDVI
            for (var i = 0; i < xCount; i++)
            {
                var index = i;
                var y1 = y;
                var offsetX = index * blockSize;
                var offsetY = y1 * blockSize;
                var actualWidth = Math.Min(blockSize, Math.Max(xSize - offsetX, 0));
                var actualHeight = Math.Min(blockSize, Math.Max(ySize - offsetY, 0));
                var bufferR = new double[actualWidth * actualHeight];
                bandR.ReadRaster(offsetX, offsetY, actualWidth, actualHeight, bufferR, actualWidth, actualHeight, 0,
                    0);
                var bufferNIR = new double[actualWidth * actualHeight];
                bandNIR.ReadRaster(offsetX, offsetY, actualWidth, actualHeight, bufferNIR, actualWidth,
                    actualHeight,
                    0, 0);
                var dataResult = new double[bufferR.Length];
                for (var j = 0; j < bufferR.Length; j++)
                {
                    var r = bufferR[j];
                    var nir = bufferNIR[j];
                    if (double.IsNaN(r) || double.IsNaN(nir)) continue;
                    var ndvi = (nir - r) / (nir + r);
                    if (double.IsNaN(ndvi)) continue;
                    var fvc = (ndvi - ndviMin) / ndviDiff;
                    dataResult[j] =
                        fvc < 0 ? NoDataValue :
                        fvc > 1 ? 1 :
                        fvc;
                }

                bandResult.WriteRaster(offsetX, offsetY, actualWidth, actualHeight, dataResult, actualWidth,
                    actualHeight, 0, 0);
            }
        }

        return (transform, xSize, ySize, srs, projection);
    }