Tuesday, September 15, 2015

Clip a raster with shapefile using C# and Gdal




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.



C# , GDAL , Raster , Vector

0 comments :

Post a Comment

 

© 2011 GIS and Remote Sensing Tools, Tips and more .. ToS | Privacy Policy | Sitemap

About Me