【C# GDAL】GDAL库、GeoTiff图像提取信息学习

1.GDAL库和GeoTIFF文件

GDAL库是处理地理信息数据的开源库。提供了一系列工具和API，供开发者能够读取、转换、写入和处理多种栅格和矢量地理空间数据格式。
GDAL提供了C++、Python、Java、C#等多种编程语言的API，使开发者能够在不同的编程环境中使用GDAL的功能。
tiff格式图片文件就是一种特殊的图片格式，与png、jpg相同。在遥感领域的影像很多都是tif格式存储的，而带有地理信息的tiff文件即可以认为是GeoTIFF

2.VS2022配置GDAL环境（C#），测试读取tiff图片

VS2022工具–NuGet包管理器–管理解决方案的NuGet程序包，直接安装GDAL包。

并且直接用应用到当前的控制台程序中。
找一张tiff格式的图片，或者用格式转换网站：https://www.zamzar.com/.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using OSGeo.GDAL;

namespace GDAL_test
{
    internal class Program
    {
        static void Main(string[] args)
        {
            //使用之前必须要配置、注册
            GdalConfiguration.ConfigureGdal();
            GdalConfiguration.ConfigureOgr();
            Gdal.AllRegister();
            string tiffFilePath = "D:\\Users\\59723\\Desktop\\su7.tiff";

            // Open函数返回了Dataset类型的对象，相当于实例化
            Dataset dataset = Gdal.Open(tiffFilePath, Access.GA_ReadOnly);

            if(dataset == null ) 
            {
                Console.WriteLine("无法打开文件");
                return; 
            }
            
            // 获取图像信息
            int rasterCount = dataset.RasterCount;
            int width = dataset.RasterXSize;
            int height = dataset.RasterYSize;
            Console.WriteLine($"宽度：{width}，高度：{height}，波段数：{rasterCount}");

            // 读取第一个波段
            Band band = dataset.GetRasterBand(1);
            if( band == null )
            {
                Console.WriteLine("无法读取波段");
                return ;
            }
            // 这个地方被坑了好久，新版本的ComputeRasterMinMax(double[] argout, int approx_ok)，已经没有out入参关键字了
            double[] minMax = { 0, 0 };
            band.ComputeRasterMinMax(minMax, 0);
            Console.WriteLine($"最小值: {minMax[0]}, 最大值: {minMax[1]}");

            // 读取波段数据
            // 在堆区开辟了一个float[]类型的数组变量，大小为width * height图片像素
            float[] rasterData = new float[width * height];
            // 参数1-2：左上角位置，参数3-4：目标区域大小，参数5：接收容器，参数6-7：容器的宽高，参数8-9：默认0
            band.ReadRaster(0, 0, width, height, rasterData, width, height, 0, 0);
            // 处理波段数据 (例如，打印前10个像素值)
            for (int i = 0; i < 10; i++)
            {
                Console.WriteLine($"像素值[{i}]: {rasterData[i]}");
            }
            dataset.Dispose();

            Console.ReadKey();
        }
    }
}

3.读取GeoTiff图片并提取相关地理信息

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using OSGeo.GDAL;

namespace GDAL_Geotiff
{
    internal class Program
    {
        static void Main(string[] args)
        {
            GdalConfiguration.ConfigureGdal();
            Gdal.AllRegister();

            string geotiff_file_path = "D:\\Users\\59723\\Desktop\\landscan-global-2022-colorized.tif";
            Dataset dataset = Gdal.Open(geotiff_file_path, Access.GA_ReadOnly);
            if(dataset == null ) 
            {
                Console.WriteLine("无法打开！");
                return;
            }

            Console.WriteLine("投影信息：");
            Console.WriteLine(dataset.GetProjection());
            //Console.ReadKey();

            double[] geoTransform = new double[6];
            dataset.GetGeoTransform(geoTransform);
            Console.WriteLine("地理变换参数：");
            Console.WriteLine("左上角X：" + geoTransform[0]);
            Console.WriteLine("像素宽度： " + geoTransform[1]);
            Console.WriteLine("旋转参数： " + geoTransform[2]);
            Console.WriteLine("左上角Y： " + geoTransform[3]);
            Console.WriteLine("旋转参数： " + geoTransform[4]);
            Console.WriteLine("像素高度： " + geoTransform[5]);

            int bandCount = dataset.RasterCount;
            Console.WriteLine($"波段数：{bandCount}");
            // 获取波段数
            for (int i=1; i<=bandCount;i++)
            {
                Band band = dataset.GetRasterBand(i); // 只有四个波段
                if (band == null)
                {
                    Console.WriteLine("无法获取波段信息");
                    return;
                }
                double min = 0, max = 0, mean = 0, stddev = 0;
                band.GetStatistics(0, 1, out min, out max, out mean, out stddev);
                Console.WriteLine($"--------波段{i}统计信息：");
                Console.WriteLine("最小值： " + min);
                Console.WriteLine("最大值： " + max);
                Console.WriteLine("均值： " + mean);
                Console.WriteLine("标准差： " + stddev);
            }
            Console.ReadKey();

            dataset.Dispose();
        }
    }
}

解释：
1.dataset.GetGeoTransform(geoTransform);
返回一个包含6个参数的数组：

• geoTransform[0]：左上角像素的X坐标（地理坐标）。
• geoTransform[1]：像素宽度（水平分辨率，表示一个像素对应的实际距离）。
• geoTransform[2]：旋转参数（通常为0，如果图像是正北向上）。
• geoTransform[3]：左上角像素的Y坐标（地理坐标）。
• geoTransform[4]：旋转参数（通常为0，如果图像是正北向上）。
• geoTransform[5]：像素高度（垂直分辨率，通常为负值，因为图像通常从左上角开始）。

地理变换参数的公式：假设你有一个像素位置 (i, j)，要计算其地理坐标 (X, Y)，可以使用以下公式：
double x = geoTransform[0] + i * geoTransform[1] + j * geoTransform[2];
double y = geoTransform[3] + i * geoTransform[4] + j * geoTransform[5];

2.band.GetStatistics(0, 1, out min, out max, out mean, out stddev);

public void GetStatistics(
    int approxOK,
    int force,
    out double min,
    out double max,
    out double mean,
    out double stddev
);

• approxOK: 是否允许使用近似值。如果为1，表示允许返回近似统计数据；如果为0，则必须计算精确值。
• force: 是否强制重新计算统计数据。如果为1，表示强制重新计算；如果为0，如果统计数据已存在，则使用现有数据。
• min, max, mean, stddev: 输出参数，分别表示最小值、最大值、均值和标准差。

4.图片元数据

元数据：为描述数据的数据，主要是描述数据属性的信息。用来支持如指示存储位置、历史数据、资源查找、文件记录等功能。
图片元数据（Metadata） 是嵌入到图片文件中的一些标签。

• EXIF：通常被数码相机在拍摄照片时自动添加，比如相机型号、镜头、曝光、图片尺寸等信息。
• IPTF：比如图片标题、关键字、说明、作者、版权等信息。主要是由人工在后期通过软件写入的数据。
• XMP：XMP实际上是一种元数据存储和管理的标准，可以将Exif，IPTC或其他的数据都按XMP统一的格式存放在图像文件中。

5.读取DEM数字高程模型的TIFF文件

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using BitMiracle.LibTiff.Classic;
using OSGeo.GDAL;
using OSGeo.OSR;

namespace ConsoleApp1
{
    internal class Program
    {
        static void Main(string[] args)
        {
            // 注册GDAL驱动
            GdalConfiguration.ConfigureGdal();
            Gdal.AllRegister();
            Gdal.SetConfigOption("GTIFF_SRS_SOURCE", "EPSG");

            // 打开DEM TIFF图像
            Dataset dataset = Gdal.Open("D:\\DEM.tif", Access.GA_ReadOnly);

            if (dataset != null)
            {
                // 获取图像的元数据
                string[] metadata = dataset.GetMetadata("EXIF");// 啥意思？

                // 打印元数据
                Console.WriteLine("Metadata:");
                foreach (string item in metadata)
                {
                    Console.WriteLine(item);
                }

                // 获取图像的投影信息
                string projection = dataset.GetProjection();
                Console.WriteLine("Projection:");
                Console.WriteLine(projection);

                // 获取图像的地理变换信息
                double[] geotransform = new double[6];
                dataset.GetGeoTransform(geotransform);
                Console.WriteLine("GeoTransform:");
                Console.WriteLine($"Origin X: {geotransform[0]}");
                Console.WriteLine($"Pixel Width: {geotransform[1]}");
                Console.WriteLine($"Rotation X: {geotransform[2]}");
                Console.WriteLine($"Origin Y: {geotransform[3]}");
                Console.WriteLine($"Rotation Y: {geotransform[4]}");
                Console.WriteLine($"Pixel Height: {geotransform[5]}");

                // 获取测区的宽度和高度
                int width = dataset.RasterXSize;
                int height = dataset.RasterYSize;
                Console.WriteLine("Area Size:");
                Console.WriteLine(width);
                Console.WriteLine(height);

                // 四波段α
                int bandCount = dataset.RasterCount;
                Console.WriteLine("波段数：");
                Console.WriteLine(bandCount);

                // 计算图像的地理边界框
                double minX = geotransform[0];
                double maxY = geotransform[3];
                double maxX = minX + width * geotransform[1] + height * geotransform[2];
                double minY = maxY + width * geotransform[4] + height * geotransform[5];

                // 创建空间参考对象
                SpatialReference srs = new SpatialReference(dataset.GetProjection());

                // 如果空间参考不是WGS84,则转换为WGS84
                if (srs.GetAttrValue("AUTHORITY", 1) != "4326")
                {
                    SpatialReference dstSrs = new SpatialReference(null);
                    dstSrs.ImportFromEPSG(4490);

                    // 创建坐标转换对象
                    CoordinateTransformation ct = new CoordinateTransformation(srs, dstSrs);

                    double[] x = new double[] { minX, maxX };
                    double[] y = new double[] { minY, maxY };
                    double[] z = new double[2];

                    ct.TransformPoints(2, x, y, z);

                    minX = x[0];
                    minY = y[0];
                    maxX = x[1];
                    maxY = y[1];
                }

                // 打印地理边界框(经纬度范围)
                Console.WriteLine("Geographic Bounding Box:");
                Console.WriteLine($"Min Longitude: {minX}");
                Console.WriteLine($"Min Latitude: {minY}");
                Console.WriteLine($"Max Longitude: {maxX}");
                Console.WriteLine($"Max Latitude: {maxY}");

                // 关闭数据集
                dataset.Dispose();
            }
            else
            {
                Console.WriteLine("无法打开DEM TIFF图像.");
            }
            Console.ReadKey();
        }
    }
}

5.1 解决warning

每次打印都会waring1这个警告，国内并没有很好的答复，更多的是忽略这个警告。
参考国外几个blog：
CRS 注册表定义问题 (EPSG:3035)
GTIFF_SRS_SOURCE 设置问题

根据英文提示：似乎是注册表中预先写好的和实际从 GeoTIFF 获得的投影不一致造成的警告。
可以去检查或打印：https://epsg.org/crs/wkt/id/4525

补充一下GIS的相关基础知识：
在国际上，每个坐标系统都会被分配一个 EPSG 代码，EPSG:4326 就是 WGS84 的代码。
空间坐标系对应的EPSG编号（中国）：https://blog.csdn.net/qq_41441896/article/details/104525296

Projected Coordinate Systems汇总：Projected Coordinate Systems

本实验数据的EPSG对应为4525：
https://epsg.io/4525

解决警告问题：
在GDAL库注册代码后加入：

Gdal.SetConfigOption("GTIFF_SRS_SOURCE", "EPSG");
// 或者用这个
// Gdal.SetConfigOption("GTIFF_SRS_SOURCE", "GEOKEYS");

使用 EPSG 可以确保投影参数与官方 EPSG 注册表一致，而使用 GEOKEYS 则会保留 GeoTIFF 文件中自定义的投影参数。

5.2 projection参数

Projection:
PROJCS["CGCS2000 / 3-degree Gauss-Kruger zone 37",GEOGCS["China Geodetic Coordinate System 2000",DATUM["China_2000",SPHEROID["CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","1043"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4490"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",1],PARAMETER["false_easting",37500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","4525"]]

• PROJCS（投影坐标系）：描述所使用的投影坐标系。
• GEOGCS（地理坐标系）：描述所使用的地理坐标系。
• DATUM（基准）：描述地理坐标系所依据的基准。
• SPHEROID（椭球体）：描述地理坐标系所依据的椭球体参数。
• PRIMEM（本初子午线）：描述本初子午线的经度。
• UNIT（单位）：描述投影和地理坐标系的单位。

PROJCS[“CGCS2000 / 3-degree Gauss-Kruger zone 37”, … ]：表示使用的是中国大地坐标系 2000（CGCS2000）的 3 度高斯克吕格第 37 带投影。

GEOGCS[“China Geodetic Coordinate System 2000”, … ]：使用的是中国大地坐标系 2000。

DATUM[“China_2000”, … ]：使用的是中国 2000 基准。

SPHEROID[“CGCS2000”,6378137,298.257222101,AUTHORITY[“EPSG”,“1024”]]：
名称：CGCS2000。
长半轴：6378137 米。
扁率：298.257222101。
AUTHORITY[“EPSG”,“1024”]：EPSG 编码为 1024。

PRIMEM[“Greenwich”,0,AUTHORITY[“EPSG”,“8901”]]：使用格林尼治本初子午线，经度为 0 度，EPSG 编码为 8901。

UNIT[“degree”,0.0174532925199433,AUTHORITY[“EPSG”,“9122”]]：地理坐标的单位是度，每度等于 0.0174532925199433 弧度，EPSG 编码为 9122。

AUTHORITY[“EPSG”,“4490”]：授权机构，EPSG 编码为 4490 的地理坐标系统（CGCS2000）。

PROJECTION[“Transverse_Mercator”]：投影类型，使用的是横轴墨卡托投影（Transverse Mercator）。

PARAMETER：投影参数
PARAMETER[“latitude_of_origin”,0]：原点纬度为 0 度。
PARAMETER[“central_meridian”,111]：中央经线为 111 度。
PARAMETER[“scale_factor”,1]：比例因子为 1。
PARAMETER[“false_easting”,37500000]：伪东移 37500000 米。
PARAMETER[“false_northing”,0]：伪北移 0 米。

UNIT[“metre”,1,AUTHORITY[“EPSG”,“9001”]]：投影坐标的单位是米，每单位等于 1 米，EPSG 编码为 9001。

AXIS[“X”,EAST]：X 轴（横坐标）指向东。
AXIS[“Y”,NORTH]：Y 轴（纵坐标）指向北。

AUTHORITY[“EPSG”,“4525”]：EPSG 编码为 4525 的投影坐标系统（CGCS2000 / 3-degree Gauss-Kruger zone 37）。

6.根据高程提取DEM边界

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using BitMiracle.LibTiff.Classic;
using OSGeo.GDAL;
using OSGeo.OSR;

namespace ConsoleApp1
{
    internal class Program
    {
        static void Main(string[] args)
        {
            // 注册GDAL驱动
            GdalConfiguration.ConfigureGdal();
            Gdal.AllRegister();
            Gdal.SetConfigOption("GTIFF_SRS_SOURCE", "EPSG");

            // 打开DEM TIFF图像
            Dataset dataset = Gdal.Open("D:\\Users\\Desktop\\DEM.tif", Access.GA_ReadOnly);

            if (dataset != null)
            {
                // 获取波段数量
                int bandCount = dataset.RasterCount;
                Console.WriteLine("波段数量: " + bandCount);

                // DEM只有一个波段存放高程数据
                Band demBand = dataset.GetRasterBand(1);
                int xSize = demBand.XSize;// 图像列宽
                int ySize = demBand.YSize;// 图像行高
                float[] demData = new float[xSize * ySize];// 建立容器存放高程数据，每个像素都有对应的高程信息
                // 读取波段信息函数ReadRaster()
                // 参数1-2：左上角位置，参数3-4：目标区域大小，参数5：接收容器，参数6-7：容器的宽高，参数8-9：默认0，不进行下采样，逐像素读取
                demBand.ReadRaster(0, 0, xSize, ySize, demData, xSize, ySize, 0, 0);

                // 创建地形掩码
                byte[] terrainMask = new byte[xSize * ySize];// byte表示八位无符号整数（1字节）
                for (int i = 0; i < demData.Length; i++)
                {
                    // 二值化处理，三目运算符处理
                    terrainMask[i] = demData[i] > 0 ? (byte)1 : (byte)0;
                }

                // 进行形态学操作以提取边界（简单实现，实际复杂应用中可以使用更复杂的算法）
                byte[] boundary = new byte[xSize * ySize];
                // 注意：该处的遍历x、y不能从0开始
                // 为了避免处理图像边界的像素
                // 如果直接从0开始，后续检查上下左右时会触发数组越界异常
                for (int y = 1; y < ySize - 1; y++)
                {
                    for (int x = 1; x < xSize - 1; x++)
                    {
                        int idx = y * xSize + x;// 定义索引，遍历像素id（第y行的第x个）
                        // 首先被检查像素一定要有高程信息(为1)
                        // 相邻上下左右四个像素，有一个没有高程信息（为0），则视为被检查像素视为边界像素
                        if (terrainMask[idx] == 1 &&
                            (terrainMask[idx - 1] == 0 || terrainMask[idx + 1] == 0 ||
                             terrainMask[idx - xSize] == 0 || terrainMask[idx + xSize] == 0))
                        {
                            boundary[idx] = 1;
                        }
                    }
                }
                // 保存结果为新的TIFF文件
                string outputPath = "D:\\Users\\Desktop\\DEM_Demo\\save_boundary.tif";
                // 获取GDAL中支持GTiff格式的驱动（用于处理TIFF文件）
                Driver driver = Gdal.GetDriverByName("GTiff");
                // outputPath: 输出文件的路径；
                // xSize: 图像的宽度ySize; ySize: 图像的高度;
                // 1: 图像的波段数，这里只有一个波段;
                // DataType.GDT_Byte: 数据类型，这里为字节类型（0-255）
                // null: 额外的创建选项，这里没有使用任何额外选项。
                Dataset outDataset = driver.Create(outputPath, xSize, ySize, 1, DataType.GDT_Byte, null);

                // 获取并设置地理变换信息,将原先DEM.tif中的数据拷贝一份
                double[] geoTransform = new double[6];
                dataset.GetGeoTransform(geoTransform);
                outDataset.SetGeoTransform(geoTransform);
                outDataset.SetProjection(dataset.GetProjection());

                Band outBand = outDataset.GetRasterBand(1);
                // 向新的DEM中写入高程信息（1和0）
                outBand.WriteRaster(0, 0, xSize, ySize, boundary, xSize, ySize, 0, 0);
                // 这行代码是去除非边界点，整个图像只显示boundary
                //outBand.SetNoDataValue(0);
                outBand.FlushCache();

                dataset.Dispose();
                outDataset.Dispose();

                Console.WriteLine("地形边界提取完成并保存至：" + outputPath);

            }
            else
            {
                Console.WriteLine("无法打开DEM TIFF图像.");
            }
            Console.ReadKey();
        }
    }
}

警告：

Warning 1: PROJ: proj_create_from_database: Cannot find proj.db

换了一台电脑后出现的问题，怀疑是没有安装好环境。暂时未解决，设计到投影信息的，对提取似乎不影响。

7.将边界换算成东、北坐标并写入Excel

using OfficeOpenXml;
using OSGeo.GDAL;
using System;
using System.Collections.Generic;
using System.Runtime.InteropServices;

namespace ConsoleApp1
{
    internal class Program
    {
        static void Main(string[] args)
        {
            GdalConfiguration.ConfigureGdal();
            Gdal.AllRegister();
            Gdal.SetConfigOption("GTIFF_SRS_SOURCE", "EPSG");

            string boundaryFilePath = "D:\\Users\\59723\\Desktop\\DEM\\save_only_boundary.tif";
            Dataset dataset = Gdal.Open(boundaryFilePath, Access.GA_ReadOnly);
            if(dataset == null )
            {
                Console.WriteLine("无法打开tif文件");
            }

            Band band = dataset.GetRasterBand(1);
            int Xwidth = band.XSize;
            int Yheigtht = band.YSize;

            // 获取地理变换参数
            double[] geoTransform = new double[6];
            // 这样取得的地理坐标数据是整个图像的左上角
            dataset.GetGeoTransform(geoTransform);

            float[] band1_data = new float[Xwidth * Yheigtht];
            band.ReadRaster(0, 0, Xwidth, Yheigtht, band1_data, Xwidth, Yheigtht, 0, 0);

            List<(double x, double y)> coordinates = new List<(double x, double y)>();
            for(int y= 0; y < Yheigtht; y++)
            {
                for(int x= 0; x < Xwidth; x++)
                {
                    int index = y * Xwidth + x;
                    if (Math.Abs(band1_data[index] - 1.0) < 0.01)
                    {
                        double geoX = geoTransform[0] + x * geoTransform[1] + y * geoTransform[2];
                        double geoY = geoTransform[3] + x * geoTransform[4] + y * geoTransform[5];
                        coordinates.Add((geoX, geoY));
                    }
                }
            }

            // 保存到Excel
            ExcelPackage.LicenseContext = LicenseContext.NonCommercial;
            using(ExcelPackage excelPackage = new ExcelPackage())
            {
                ExcelWorksheet worksheet = excelPackage.Workbook.Worksheets.Add("sheet1");
                worksheet.Cells[1, 1].Value = "东坐标";
                worksheet.Cells[1, 2].Value = "北坐标";
                for(int i= 0; i < coordinates.Count; i++)
                {
                    // 设置单元格格式，保留两位小数
                    worksheet.Cells[i+2, 1].Style.Numberformat.Format = "0.000000000";
                    worksheet.Cells[i+2, 2].Style.Numberformat.Format = "0.000000000";
                    worksheet.Cells[i+2, 1].Value = coordinates[i].x;
                    worksheet.Cells[i+2, 2].Value = coordinates[i].y;
                }
                string save_file = "D:\\Users\\59723\\Desktop\\DEM\\Coordinates.xlsx";
                excelPackage.SaveAs(save_file);
                Console.WriteLine("文件保存至：" + save_file);
            }
            Console.ReadKey();
        }
    }
}

将边界坐标信息成功保存到Excel文件中，使用图标粗略展示，和arcgis中展示的相同。