当前位置：首页 > news >正文

GDAL+C#实现矢量多边形转栅格

news 2024/10/22 14:41:06

1. 开发环境测试

参考C#配置GDAL环境，确保GDAL能使用，步骤简述如下：

创建.NET Framework 4.7.2的控制台应用

注意：
项目路径中不要有中文，否则可能报错：can not find proj.db

在NuGet中安装GDAL 3.9.1和GDAL.Native 3.9.1
测试代码（读取一张影像的长、宽、波段数、坐标系）

using OSGeo.GDAL;
using System;namespace TestGDAL
{internal class Program{static void Main(string[] args){Configure();Test();Console.ReadLine();}public static void Configure(){GdalConfiguration.ConfigureGdal();GdalConfiguration.ConfigureOgr();Gdal.AllRegister();}public static void Test(){Dataset ds = Gdal.Open(@"./测试.tif", Access.GA_ReadOnly);int rasterX = ds.RasterXSize;//影像宽度int rasterY = ds.RasterYSize;//影像高度int bandCount = ds.RasterCount;//波段数double[] tmpD = new double[6];ds.GetGeoTransform(tmpD); //影像坐标变换参数//影像坐标系信息（WKT格式字符串）string proj = ds.GetProjection();//用一个消息弹窗把信息显示出来Console.WriteLine($"宽度={rasterX}，高度={rasterY}，波段数={bandCount}");Console.WriteLine($"坐标系={proj}");}}
}

2. 多边形转栅格（按包围盒区域）

参考：C#调用GDAL实现矢量转栅格

用ArcMap创建一个shp面数据，用于后续输入。
基于投影坐标系WGS_1984_UTM_Zone_50N，绘制了三个面，并添加属性字段Code，该字段将作为栅格数据的值。
示例代码

using OSGeo.GDAL;
using OSGeo.OGR;
using System;
using OSGeo.OSR;
using Driver = OSGeo.GDAL.Driver;namespace TestGDAL
{public class Test2{static string shpPath = @"./polygon.shp";public static void Run(){//初始化GDAL和OGRGdalConfiguration.ConfigureGdal();Gdal.AllRegister();Ogr.RegisterAll();// 获取矢量图层DataSource vectorDS = Ogr.Open(shpPath, 0);Layer vectorLayer = vectorDS.GetLayerByIndex(0);//设置栅格化参数，按1000行1000列来栅格化Envelope extent = new Envelope();vectorLayer.GetExtent(extent, 0);double xMin = extent.MinX;double xMax = extent.MaxX;double yMin = extent.MinY;double yMax = extent.MaxY;int xSize = 1000;int ySize = 1000;double xRes = (xMax - xMin) / xSize;double yRes = (yMax - yMin) / ySize;Console.WriteLine("Extent: " + extent.MaxX + " " + extent.MinX + " " + extent.MaxY + " " + extent.MinY);Console.WriteLine("X resolution: " + xRes);Console.WriteLine("Y resolution: " + yRes);// 创建栅格化输出文件var time = DateTime.Now.ToString("yy-MM-dd HH-mm-ss");string rasterFile = $"./output[{time}].tif";Driver rasterDriver = Gdal.GetDriverByName("GTiff");Dataset rasterDS = rasterDriver.Create(rasterFile, xSize, ySize, 1, DataType.GDT_Float32, null);// 设置栅格化文件参考系SpatialReference rasterSRS = vectorLayer.GetSpatialRef();rasterSRS.ExportToWkt(out string wkt, null);rasterDS.SetProjection(wkt);rasterDS.SetGeoTransform(new double[] { xMin, xRes, 0, yMax, 0, -yRes });Console.WriteLine("Projection:\n " + wkt);//设置属性字段为Codestring[] rasterizeOptions = new string[] { "ATTRIBUTE=" + "Code" };//执行栅格化Gdal.RasterizeLayer(rasterDS, 1, new int[] { 1 }, vectorLayer, IntPtr.Zero, IntPtr.Zero, 1, new double[] { 1 }, rasterizeOptions, null, null);//将结果输出并关闭文件rasterDS.FlushCache();rasterDS.Dispose();vectorDS.Dispose();}}
}

栅格化效果图

在这里插入图片描述

3. 多边形转栅格（只保留多边形区域内的栅格）

在前面基础上，如果只想保留多边形内的栅格，也就是想要下面的效果：

在这里插入图片描述
基于上面的测试数据，则可以通过调整栅格化参数来实现：

//执行栅格化
rasterDS.GetRasterBand(1).SetNoDataValue(0);//0值表示没有值
Gdal.RasterizeLayer(rasterDS, 1, new int[] { 1 }, vectorLayer, IntPtr.Zero, IntPtr.Zero, 0, new double[] { 0 }, rasterizeOptions, null, null);