C# 讀取四波段遙感影像生成植被覆蓋度柵格(TIF)

turingguo發表於2024-11-15
C#

GDAL使用

使用gdal.netcore來讀取和生成柵格檔案。
優點：自帶gdal執行時相關檔案，不用額外再安裝gdal庫
缺點：導致釋出的檔案變大很多，比如Win+Linux的執行時加起來就超過了400M，所以最好是按需載入對應的執行時。比如DEBUG是執行在win的，就新增MaxRev.Gdal.WindowsRuntime.Minimal包，然後Release環境是執行在linux x64的，就新增MaxRev.Gdal.LinuxRuntime.Minimal.x64。

系統 執行時
windows MaxRev.Gdal.WindowsRuntime.Minimal
linux-x64 MaxRev.Gdal.WindowsRuntime.Minimal
linux-arm64 MaxRev.Gdal.WindowsRuntime.Minimal
mac-x64 MaxRev.Gdal.WindowsRuntime.Minimal
mac-arm64 MaxRev.Gdal.WindowsRuntime.Minimal

如圖所示，僅支援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);
    }

相關文章