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);
}