How to clip raster data to a polygon using GDAL and OGR in C#

Rare work-related post: Map data manipulation using older GDAL and OGR libraries in C#, which is a more obscure API compared to many Python and CLI answers found online.

NuGet packages: GDAL 1.11.1, GDAL.Native 1.11.1

using OSGeo.GDAL;
using OSGeo.OGR;

// Provided parameters
// For clarity, these are declared explicitly
// But they should be passed as arguments to a method

// Path to original raster image, shape, and rasterized polygon image
var rasterFilePath = "C:/Temp/image.tiff";
var polygonFilePath = "C:/Temp/polygon.shp";
var rasterizedPolygonFilePath = "C:/Temp/polygon.tiff";
// Extent / bounding box for the image
// As an example, a field in Ohio
var north = 41.60329;
var south = 41.58917;
var east = -83.20192;
var west = -83.18218;
// Which value on the image means "no data"
// As an example, solid white
var noDataValue = 1.0;
// Polygon in OGR format
// As an example, created from some unspecified JSON
OSGeo.Geometry.Polygon clippingPolygon = Ogr.CreateGeometryFromJson("{}");

GdalConfiguration.ConfigureGdal();
Ogr.RegisterAll();

// Open raster image dataset to be updated
using(var rasterDataset = Gdal.Open(rasterFilePath, Access.GA_Update))
{
	// Rasterize the clipping polygon
	using (var polygonDataSource = Ogr.GetDriverByName("ESRI Shapefile").CreateDataSource(polygonFilePath))
	{
		// Create an OGR layer and add the polygon as a feature to it
		var polygonLayer = polygonDataSource.CreateLayer("Polygon", null, wkbGeometryType.wkbPolygon, null);
		var polygonFeature = new Feature(polygonLayer.GetLayerDefn());
		polygonFeature.SetGeometry(clippingPolygon);
		polygonLayer.CreateFeature(polygonFeature);

		// Create a new TIFF file with the same size as the original and rasterize the polygon into it
		using (var rasterizedPolygonDataset = Gdal.GetDriverByName("GTiff")
			.Create(rasterizedPolygonFilePath, rasterDataset.RasterXSize, rasterDataset.RasterYSize, 1, DataType.GDT_Byte, null))
		{
			var pixelWidth = (east - west) / rasterDataset.RasterXSize;
			var pixelHeight = (north - south) / rasterDataset.RasterYSize;

			// Your original raster image might be north-down (since screen down and north are both a positive direction)
			// Make sure to set the geotransform correctly
			
			// Per GDAL docs, GeoTransform array contains:
			// [0] X-coordinate of the upper-left corner of the upper-left pixel of the image.
			// [1] W-E pixel width
			// [2] row rotation (typically zero)
			// [3] Y-coordinate of the upper-left corner of the upper-left pixel of the image.
			// [4] column rotation (typically zero)
			// [5] N-S pixel width
			var rasterizedPolygonGeoTransform = new double[] { west, pixelWidth, 0, south, 0, pixelHeight };
			rasterizedPolygonDataset.SetGeoTransform(rasterizedPolygonGeoTransform);

			Gdal.RasterizeLayer(rasterizedPolygonDataset, 1, new int[] { 1 }, polygonLayer, IntPtr.Zero, IntPtr.Zero, 0, null, null, null, null);
		}  // Close the rasterized polygon dataset to write it to disk
	} // Close the polygon datasource, done using it

	// Re-open the rasterized polygon
	using (var rasterizedPolygonDataset = Gdal.Open(rasterizedPolygonFilePath, Access.GA_ReadOnly))
	{
		// Mask the original TIFF file by filling all points outside the polygon with No Data value
		var dataBand = rasterDataset.GetRasterBand(1);
		var maskBand = rasterizedPolygonDataset.GetRasterBand(1);

		dataBand.SetNoDataValue(noDataValue);
		Gdal.FillNoData(dataBand, maskBand, 1.0, 0, null, null, null);
	} // Close the rasterized polygon Dataset, done using it
} // Close the original image to write it to disk