C# 读取四波段遥感影像生成植被覆盖度栅格(TIF)
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);
}
