This is the second part of a blog series on mapping circular fields. In part one I talked about the challenges for mapping in Desert environments and about how thresholding and the Watershed Algorithm can be used to detect coins – this offers a potentially useful way to map circular fields. These are a challenge for planning a seismic survey whether it be through permitting or the actual operations of a seismic crew.
Satellite Data
Last week I drew on a simple cluster of coins and created a boundary around them. This time I want to draw around the circular fields in the above image. I then want to be able to access this data as a vector and use it in mapping software, primarily GIS but not exclusively as increasingly more software can read ‘geospatial’ data.
Download a Landsat8 scene or Sentinel2a scene. There are plenty of ways to do this, for ease of use I like to use the Semi Automatic Classification Plugin for QGIS. It has the added advantage that it will perform all the necessary corrections for you in post processing.
This tutorial is going to use a mix of QGIS, Python and GDAL.
Step 1.
In QGIS build a Virtual Raster Catalog (vrt). Raster – Miscellaneous – Build Virtual Raster (Catalog). Select your bands (at least 3 for an RGB). You should end up with a composite image. To speed things up I am going to zoom in on my data and clip to the extent of my area of interest (you can see this in the image above – it’s not a full scene). Clipping a raster is again easy in QGIS, Raster – Extraction – Clipper (options for Mask or extent, choose the best for you).
Step 2.
Once finished, fire up the command line with access to GDAL (OSGeo4W Shell is ideal) for this. Navigate to your newly cropped image and run the command
gdal_translate -b 1 -b 2 -b 3 -of JPEG -scale -co worldfile=yes out.tif output.jpg
Where out.tif is your clipped satellite image and output.jpg is the image to be processed. Not only are we converting to jpg but we are also writing out a world file to retain its positional information as this is needed to create the vector data. I am also just selecting 3 bands from my clipped image.
At this point we have an image that we can process in the watershed code. It will strip out the georeferencing but we can put that back at the end.
Step 3.
Run the watershed code. Here is the thresholded image
and here is the image with the boundaries mapped, they are shown in red. Looking good 🙂
Step 4.
Remember that world file created in step 2? Copy and paste it and change the name to your newly created image with the red boundaries on. This will allow you to load it into the correct place in your GIS of choice. Great! Do this for the thresholded image as well, as this is the one that is ideal to convert to a vector.
Step 5.
Now I need to get the image to a binary one so I can convert it to vector. I’ve previously written about that here
http://www.acgeospatial.co.uk/blog/reclassify-and-numpy/
You can download this code from Github. Run the reclass.py on your thresholded image.
Step 6. Final step
Back to the command prompt and GDAL. This time run the following call
gdal_polygonise reclass.tif out.shp
Where the reclass.tif image is the one from step 5 and out is the shapefile to be created.
Result
Ultimately this shapefile layer will need cleaning up a bit. The thresholding values can be changed in the Watershed Algorithm to tune it to the area you are working in. However the need to digitise is significantly reduced, attributes can be assigned to these polygons for ownership, exclusion zones created. A new circular field appeared? Get the latest image and run the process again!
Final part next week, I’ll share the code and talk about Automation.