Fun with raster clipping

Hi there and welcome to the second blog in this series where we try to determine the volume difference of excavated material in one of the biggest open pit mines, the Bingham Canyon Mine. In the first blog of this series, we have selected a really old map (thank you, open access USGS) of the area, all the way from 1896, just in the beginning of the excavation in the area. Then the map was georeferenced and visualized on a web-based platform GoldenAI GUI. If you are just joining and are interested in this step, please check out the first blog here.

And that brings us to today's topic of preparing our georeferenced raster for further analysis. Before a DEM of old times can be made and compared to a modern satellite DEM, we need to get rid of that black and white edge around the map. 

Now that's simple enough, first you go back to QGIS, then you check with the Identify tool the band values of the pixels you want to make transparent. From there, you navigate to the Transparency in Layer Properties and add the band values to the Transparency Pixel List (after going back and forth a couple of times because you keep forgetting the band values). And that's it.

Not great, not terrible, but we can make that better. All these steps were great to get rid of that black edge, which could work great if you have, for example, a satellite image, but in our case, there is a full scanned map along with all the margins. And we don't want or need that. So let's get rid of it. The same workflow is also applicable if your area of interest is smaller than the full extent of the raster and you want to reduce it for more precision and/or faster analysis. 

The simplest way to make this happen is to clip your raster. Just one quick note, before we go any further, all the examples shown are made in QGIS version 3.32., so older and, in the future, newer versions of QGIS may or may not have the same tool at the same place. 

With that disclaimer out of the way, let’s clip a raster. The tool to achieve this is in the Raster toolbar, under Extraction. There are two ways to clip a raster: Clip Raster by Extent and Clip Raster by Mask Layer. The main difference is that the first doesn't require any additional files to run, while the other, you guessed it, does. The second also gives you more freedom in choosing the shape of your output raster, but we’ll get to that later. 

In Clip Raster by Extent there are several option to what can be used as Clipping Extent:

• Calculate from Layer, Calculate from Layout Map and Calculate from Bookmark → uses the extent of another layer/layout/bookmark in the working project, great to know, but not helpful in our case

• Use Current Map Canvas Extent → zoom on a map layer so only what you want in your clipped output raster is visible; great and very simple possibility, but can be problematic if it is hard to fit what you want on a map view 
• Draw on Map Canvas → click and drag a rectangle on the input raster; again great and simple possibility, with the same potential problem if you area of interest can not neatly be placed in a rectangle

Just to get the feel for it, let's do Use Current Map Canvas Extent. After selecting it as Clipping extent, only other required input is of course the raster you want to clip. Make sure you have zoomed on the Map Canvas to the area you want to keep. There are other additional parameters that can be set, and one that I highly recommend is Assigning a Specific Nodata value to Output Bands. Why and how it is important is a bit below in the text, but for now, just know it is highly recommended to set it and to use some value like -9999, or something you are sure is outside of your value range. 

And with that, we can run the tool. If the Sun and retrograde Mercury are aligned, that is the end of it and you can skip to Figure 9. However, there is a possibility you will run into an error like this, same as I did. 

The cause of this error message is simple actually: the option Save to Temporary File is selected automatically and can work great. It produces the wanted results, adds it to your Layers and you can view it before exporting it to be sure you got what you wanted and not play “I’m feeling lucky” like in Google search. However, sometimes the file format can be automatically set to -of GMT, which does not support multiband raster clip.

Good thing there are simple solutions: one way to edit the command line and second is to just save the file as a GeoTIFF with Save to File. Both of these options end up with the GDAL/OGR console call changed to -of GTiff (which has the ability to clip multiband rasters). See, A B C it’s easy as 1 2 3 (please feel free to have the Jackson 5 stuck in your head for the rest of the day). 

After that little sidetrack, what we end up with, after running the clip tool, is something like this (the background semi transparent full map is there just to help you visualize where we are).

Now that we covered how to clip by a rectangle, let's say you want/need to clip by some other shape that more precisely represents your area of interest. That can be done in almost all the same steps with Clip Raster by Mask Layer. As mentioned in the beginning, this one requires an additional polygon layer as a mask. Sometimes it is referred to as a cookie cutter, meaning everything inside the mask (cookie cutter) is what you want (cookie), and everything outside is not important to you.

Here there are a bit more things to keep in mind. First, you need your input polygon mask layer. It can be something specific you already have or you can create a new layer and draw a polygon. Be careful that the raster (input layer) and polygon (mask layer) have the same coordinate system. There is an additional possibility to change the target coordinate reference system (in which coordinate system will your output layer be in), if you need it to be in a different CRS than the input layer. 

Second, let's go back to No Data values. You can again set the No Data value to -9999 like before, but now there is a possibility to create an output alpha band. What is an alpha band, you ask? Well it is just a binary transparency mask with a transparency value for each pixel. On the surface, both options give the same result, but under the hood, adding the alpha band does just what its name say (don’t I love simple descriptive names): transforms your 3-band raster to a 4-band raster, exactly what we need to work with OCLI. So that any further processing, like reprojection, doesn’t cause the black edge from the beginning to be introduced again. On the contrary, just specifying No Data values, without adding the alpha band, will at first glance work the same, but the dreaded black edge can be reintroduced with reprojection and some further analytical processes can be a bit harder if the alpha band is not present. So my humble recommendation is to do both, set the no data value and enable the Create an Output Alpha Band option. 

Third, make sure you also enable Match the Extent of the Clipped Raster to the Extent of the Mask Layer and the final fourth step is to Save to a GeoTIFF file.

Et voilà, you clipped the raster with the polygon of your choice. Now I know that doesn't look much different than the previous clip (either by map extent or draw on map canvas), but again, it gives you much more freedom if your area of interest is “weirdly shaped”. Sometimes it's a watershed extent, sometimes its province boundaries, sometimes you just want to play. So if the rectangle just doesn't work for you, Clip by Mask Mayer has you covered. 

Now, I didn't make this cat map just to amuse myself (maybe just a little). I actually wanted to show you what happens if you don't set a No Data value. You end up with something like this, a black rectangle around the clipped raster. It can be fixed with additional steps in the transparency layer, bringing us back to step one of this blog, but why give yourself more work if you don't have to. Speaking of giving yourself more work than necessary, if you don't have a specific irregular area of interest or a specific need to clip to anything other than a rectangle, one of the simpler clip functions without the mask is a way to go for you. 

Let’s recap what was done up until now: first we selected, georeferenced and confirmed the result in GoldenEye GUI and now we have clipped the unnecessary margins on the georeferenced map. 

And our next steps are to make a DEM out of the georeferenced and clipped old map and calculate the height and volume difference between then (old map turned to DEM) and now (new satellite image/DEM). Hope to see you soon and until then have fun with clipping your rasters. 

‍