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);}//计算NDVIfor (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%对应的Countvar 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++){//计算NDVIfor (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);}
