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.



Comments

Post a Comment

Popular posts from this blog

Prevent WPF Global Keyboard Hook from stops working after hitting keys a while C# .NET

Image
The article written by Dylan Currier,  Low Level Global Keyboard Hook / Sink in C# .NE T is one the simple, easily understandable, and and working post on Global Keyboard hook out there in internet. It worked in my WPF application and I added following line in the WPF app to refresh hook after each key press- IDK  if it is write or wrong, more likely to be a wrong practice to hook/unhook every time but it gets the job done. The primary reason to refresh hook is to prevent app hanging (not responding to the key press) after few clicks on the WPF button. void _listener_OnKeyPressed(object sender, KeyPressedArgs e)         {                       _listener.UnHookKeyboard();           _listener.HookKeyboard();         } Following is the video from  Dylan Currier   on  Low Level Keyboard Hook in C# / .NET"
Read more

Generate ArcGIS Token By URL Request From ArcGIS Portal, Federated Environment , ESRI JS API and Angular snippets

In my custom angular web application, I had to pull the data from different ArcGIS server environments, say Dev ArcGIS server and Test ArcGIS server environment, some of my items are directly referred from ArcGIS server services and some of them are from Enterprise portal referring the WebMapId and PortalId from both environment (Dev, Test). Servers are in federated environment. To pull the data from different ArcGIS server environment, user must login in each environment. In addition, user must login to get inside the custom application and Enterprise AD Group was set to authenticate users on the custom web application and ArcGIS server environments. So, there will be 3 login attempts (1 –application itself, 2- Dev server/portal, 3- Test server/portal) to use the application, which doesn’t provide good user experience. To improve the better application workflow, I decided to reduce the number of logins required in the application and use AD Username and Password captured du...
Read more

Solution: IntelliJ not recognizing a particular file correctly, instead its stuck as a text file

Image
I had mistakenly save my class as fileName.text file and tried to rename it with fileName .java but the IntelliJ couldn't recognize it and throws the error saying.. After spending long hours on internet, I found a solution and sharing/keeping here for future reference. Step 1: Click "File"==> "Settings" Step 2: Expand "Editor" & Click "File Types" Step 3: You will see all file types on Right. Navigate to the "Text Files" and Click it Step 4: You should able to see your file name on the bottom of Registered Patterns (lower box) Step 5: Remove your file from the Registered Patterns. The problem should solved and let you to rename with fileName.java Step 6: If not, delete the file from the project and create it again with name fileName
Read more