Over the week, I stumbled on problem clipping raster with feature using C# and Gdal and come up with following solution, which clips the features and fills no data value for missing raster values inside the extent. If you are creating a new project- set up GDAL & C# environment as described here
1: string featureName = "myshp.shp";
2: string rasterName = "myraster.tif";
3: int rasterCellSize = 10;
4:
5: //Register vector drivers
6: Ogr.RegisterAll();
7:
8: //Reading the vector data
9: DataSource dataSource = Ogr.Open(featureName, 0);
10: Layer layer = dataSource.GetLayerByIndex(0);
11:
12: //Extrac srs from input feature
13: string inputShapeSrs;
14: SpatialReference spatialRefrence = layer.GetSpatialRef();
15: spatialRefrence.ExportToWkt(out inputShapeSrs);
16:
17: Envelope envelope = new Envelope();
18: layer.GetExtent(envelope, 0);
19:
20: //Compute the out raster cell resolutions
21: int x_res = Convert.ToInt32((envelope.MaxX - envelope.MinX) / rasterCellSize);
22: int y_res = Convert.ToInt32((envelope.MaxY - envelope.MinY) / rasterCellSize);
23:
24: //Register vector drivers
25: Gdal.AllRegister();
26:
27: Dataset oldRasterDataset = Gdal.Open(rasterName, Access.GA_ReadOnly);
28:
29: //Create new tiff
30: OSGeo.GDAL.Driver outputDriver = Gdal.GetDriverByName("GTiff");
31:
32: //New geotiff name
33: string outputRasterFile = "mynewraster.tif";
34:
35: Dataset outputDataset = outputDriver.Create(outputRasterFile, x_res,y_res, 1, DataType.GDT_Float64, null);
36:
37: //Geotransform
38: double[] argin = new double[] { envelope.MinX, rasterCellSize, 0, envelope.MaxY, 0, -rasterCellSize };
39: outputDataset.SetGeoTransform(argin);
40:
41:
42: //Set no data
43: Band band = outputDataset.GetRasterBand(1);
44: band.SetNoDataValue(-9999);
45: outputDataset.SetProjection(inputShapeSrs);
46:
47: string[] reprojectOptions = {"NUM_THREADS = ALL_CPUS"," INIT_DEST = NO_DATA","WRITE_FLUSH = YES" };
48:
49: Gdal.ReprojectImage(oldRasterDataset, outputDataset, null, inputShapeSrs, ResampleAlg.GRA_NearestNeighbour, 1.0,1.0, null, null, reprojectOptions);
50:
51: //flush cache
52: outputDataset.FlushCache();
53: band.FlushCache();
54: dataSource.FlushCache();
I am replacing ESRI ArcGIS dependencies with Gdal as much as I can for raster processing. If anybody has elegant solution please feel free to suggest.

0 comments :
Post a Comment