My blog, imported from Blogger and converted using Jekyll.
I have been wondering about how it might be possible to automatically decide on land cover classification rules, as I mentioned in previous posts.
I realised the first step was to create a way of using a list of the images of the ground reference points and their coordinates, to examine the images together with remote-sensed images.
I've made a Python program that takes in a list of image file names and their latitude and longitude co-ordinates from a csv file, and uses these to make a subset of an area of a Landsat image 3km x 3km, and display a spectrum of a smaller area around the point. I presently use an area of 120m x 120m which is 4x4 Landsat pixels.
The code expects a layerstacked image of the optical and near-infrared bands, of which there are 6 for Landsat 7, and 7 for Landsat 8.
It will convert the lat/long coordinates, to OSGB grid references for labelling, and if specified to projected UTM coordinates which the Landsat images are released in.
The output below uses a Landsat 8 image in UTM30N, it is labelled in OSGB, I believe UTM30N grid north is fairly close to OSGB36.
Starting with a few of the pictures from the MSc fieldtrips I mentioned:
The images below are from Landsat 8 in July 2014, with the fieldtrips in Feb/March 2014. The versions in the Google+ album below are from Landsat 7 data from March 2007.
There is a little caveat I will make that certain of the ground level pictures are looking at a point a distance away from the coordinates that are where the photographer was standing.
Also the GPS sometimes made catastrophic errors, and occasionally I have manually edited the coordinates by taking an average of other pictures at the same site (within a small area).
The code is now available on my Bitbucket account
. It only works on Linux and is likely to still be a bit rough and ready at the time of writing.
I have also uploaded the output of the fieldwork ground reference points as a Google+ album
|On the Penglais campus at the back entrance of the Llandinam building.|
|Ynyslas dunes near Borth|
|An area of decidious woodland near Tal y Bont|
|Some upland grassland with Dr Pete Bunting|
|Vaccinium near Nant y Arian|
|Looking west, towards Taylor's Leat, an ingenious method of providing water to drive a wheel to power mining machinery in the 19th century.|
|Next to Nant y Moch reservoir dam|
The Department of Work and Pensions insists that all jobseekers use their online jobsearch website, 'Universal Jobmatch
Unfortunately, it can be quite difficult to use, with needing to click through a load of pages filled with mostly agency adverts, which are often duplicates, and sometimes cross-posted to incorrect locations despite the website offering an option to search within a particular radius.
So I thought I'd write a Python script to convert the information to a simple HTML table that one can simply scroll down through without needing to click through 20 or so pages.
I have now tidied it up a bit at put it on Bitbucket: Universal Jobmatch Spam Soup
. As I used the Beautiful Soup
library and it is written in Python I thought 'Spam Soup' was an appropriate name.
It doesn't really do much to the results, in the way of filtering them beyond what the query returns, but it does make it easier to read I think.
In the segmentation discussed in the previous two posts, as well as mean elevation, slope and aspect, I also have two curvature layers, longitudinal, and cross-sectional curvature that I calculated in LandSerf
. The definitions of these are explained in the PhD thesis of Jo Wood, the author of LandSerf
I thought to indicate longitudinal curvature using colour, shading convex slopes in red and concave in blue. Cross-sectional curvature is indicated with dots, which are white for diverging slopes and black for converging. Both of these are normalised by the square root of the segment area relative to the minimum segment area.
|Newquay in Cornwall. Convex slopes appear in red, concave in blue.|
|Cadair Idris. Some strongly curved areas create very large dots that fill entire segments above Llyn Cau|
|Rescaling the dots to half their radius.|
|Normalising the radius of the dots by the slope. i.e. reduce dots on steep slopes, and expand those on flatter terrain, such that they are scaled by the magnitude of the cross-sectional curvature relative to the slope. It does get a bit complicated, with the formula for their size in QGIS being "case when RAT_CrScCrv_Mean > 0 then 5*(90/max(RAT_slope_deg_Mean,1))*RAT_CrScCrv_Mean*(sqrt($area/(14*24*24))) else -5*(90/max(RAT_slope_deg_Mean,1))*RAT_CrScCrv_Mean*(sqrt($area/(14*24*24))) end"|
This is another map like the ones in yesterday's post
, this time of the town of Aberystwyth in Wales, and a little of the area surrounding it.
I have altered the spacing of the slope lines, which are now 20m apart, with grid squares shown at 1km size.
Roads as usual come from OpenStreetMap, with OS Vector Map foreshore, woodland, plus tidal and inland water (the labels of these are from OpenStreetMap - possibly a dangerous thing to do, if something is in OSM but not OS VectorMap there will be a labelled object that is not mapped!).
Hill peaks are from www.hill-bagging.co.uk
, plotting all greater than 150m with more than 20m drop.
As before, the steepness of hills is indicated by the width of the slope lines, in this version, they are plotted with 10*((Slope_deg_Mean - 2.25) / 45.0) metres width, if the slope exceeds 2.25 degrees, and if not omitted.
That is no lines are plotted where the average slope of the segment is less than 1 in 20, and increase in width up to 19.5m as slope approaches 90 degrees. So a vertical cliff will be almost solid red.
They are plotted in the direction of the mean aspect of the segment. This is the real aspect, not the absolute distance from north.
It is possible that a segment facing north may for example have half of it at 1 degree east, and the other half at 359 degrees. Thus the average would be at 180 degrees. This doesn't matter because the lines have no arrows on them.
However if you have 3/4 of it facing 359 degrees, and 1/4 facing 1 degree, then it would produce spurious results, therefore slope direction may be dubious on north-facing slopes.
|Aberystwyth, with the Ystwyth and Rheidol rivers, and Clarach Bay at the top. There are certain areas where a north-facing slope is given a spurious east-west slope line, such as the north-facing slope at the top of this image at Clarach.|
When I originally did this, I thought that the excess division of segments along the N/S line when I segmented using the original aspect, rather than the absolute angle from north, was a problem. However for this purpose, dividing segments that cross the 360/0 degree discontinuity helps some of the north-facing slopes avoid spurious east-west 'mean aspect' of segments:
|Compare the north-facing slopes to the previous image.|
|The segment boundaries plotted for illustrative purposes.|
Here I demonstrate a scheme of visualising slopes on a map using lines to indicate direction and steepness of slope.
This is not necessarily easy, because both vary continuously across a terrain.
So therefore I thought to use RSGISLib
to segment a digital elevation model, that I had derived from SRTM 1 arcsec data, and used GDAL
to create derived topographic variables and make a layerstack that consisted of 6 bands:
- Elevation above sea level (m)
- Slope (degrees)
- Aspect (degrees)
- magnitude of Aspect away from N (degrees) - to avoid the 360/0 discontinuity
- Cross-sectional curvature
- Longitudinal curvature
I told RSGISLib to use all bands except band 3 - to avoid the 360/0 discontinuity dominating the segmentation and getting a large number of segments divided along the N/S axis.
After producing the segmented DEM, then it is a matter of using RSGISLib to populate the Mean and StdDev of each band to the raster attribute table of the segment clumps file, and creating a shapefile version using gdal_polygonize and using RSGISLib to export the raster attribute table into a CSV file. There is a fair bit of scripting involved, but I was able to recycle the scripts from my MSc dissertation.
Then I load this up into QGIS, and create some styling rules to show the lines in the direction of the aspect, and thickness according to the steepness of the slope.
I added a colourisation of the elevation to display that as well. I have neglected to show a key for this, but qualitatively the colour ramp begins with a pale green becoming stronger with increasing altitude up to 200m, with yellow overlaid gradually coming in above 100m up to 300m, then orange/brown comes in up to the maximum of 420m.
Here are a couple of areas of Cornwall mapped according to this scheme. I have shown segment boundaries lightly (probably only visible on the png versions linked to in the captions) which have a minimum size of 14 pixels, where the pixel size is approximately 24m.
|A larger version, showing a wider area including Redruth/Camborne (14MB PNG).|
Part of N coast and Bodmin Moor
Close up of Truro area, rendered at 1:25000
Here is a version, with the roads added from OpenStreetMap, and inland and tidal water added from OS VectorMap, rendered at 1:25k with a modified spacing of slope lines which now are removed entirely where the mean slope in a segment is less than 2.25 degrees (1 in 20), and plotted in a shade of red in an effort to avoid confusion with minor roads. The segment boundaries are also no longer plotted.
|A larger version, rendered at a higher resolution, covering a somewhat wider area (7.6MB PNG).|
Redruth, and the north face of Carn Brea
|The elevation scale is modified such that darker green is lower in this rendering, moving to light green, through yellow and then orange as before. The north face of Carn Brea shows areas of steep slope, which can sometimes suffer from gorse fires. |