Natural Resources Database - GIS Files

Raster Analyser v0.5

Written by Richard Alexander, 29th June 2012

Purpose

The Raster Analyser tool enables you to rapidly calculate area intersects between two or more GIS files using raster processing.

For example you can calculate the:

• Total area of polygons from one layer intersecting individual polygons in a second layer.

• Combined area of polygons (with overlaps removed) from multiple layers intersecting polygons in a second layer.

Features:

• Input files can be in ESRI shapefile format or MapInfo table format.

• Supports multiple input or intersect layers.

• Input layers may be combined, intersected or processed individually.

• Intersect layer polygons may be grouped by attribute values

• Raster processing is performed on a grid of 10m x 10m pixels by default.

• Intersect areas are calculated in hectares and output to a text file

Using the tool

Run the tool by opening the RasterIntersect.exe file.

Input Layer

Add the layers which are to be intersected and their areas calculated. These can be ESRI shapefiles, MapInfo tables or MapInfo MID/MIF files.

If Combine input layers is selected then the layers will be merged and the footprint of the combined layers will be processed.

If Intersect Input Layers is selected then only the intersecting area of the input layers will be processed.

If _Separately process input layers _is selected then the results will be shown for each input layer.

Intersect layers

Add the layers that will be used to intersect the input layers. These can be shapefiles or MapInfo files.

If more than one layer is added then these will be merged into a single layer that will be used to intersect the input layers.

If only one layer is added you can choose to group the results by selecting one of the attributes.

If grouping you can choose to add all the intersect areas for an attribute by selecting Sum, or remove duplicate rows for that attribute, by selecting First.

Start processing

When you click on ‘Start processing’ each of the input layers will be loaded in turn and then intersected with the intersect layers. If a Group by attribute has been selected then the input layers will be intersected against the individual polygons rather than the entire layer.

A table will open showing the results:

The table is in tab separated format that can be opened in a Excel.

Example query:

A) Calculating the extent of lowland grassland within Sites of Scientific Interest/Higher Level Stewardship

i) Add the 5 lowland grassland inventories (lowland calcareous grassland, lowland dry-acid grassland, lowland meadows, purple moor grass and rush pasture, upland hay meadow) to the Input layer.

Select Combine input layers from the drop list

ii) Add the HLS and SSSI layers to Intersect layers.

iii) Click on Start Processing…

The resulting table will show the combined area of the lowland grassland inventories and the area of this that occurs within SSSIs/HLS land (i.e. within either SSSI or HLS land or within both).

If you want a breakdown of the area of lowland grassland exclusively within SSSIs, exclusively within HLS and within both, re-run the query with the SSSI layer on its own and then the HLS layer on its own. You can then calculate the breakdown of SSSI/HLS/both from the three results.

B) Calculating the area of each priority habitat for each SSSI

i) Add each priority habitat inventory to the Input layers.

Make sure that Process each layer separately is selected.

ii) Add the SSSI layer to the Intersect layers.

Select the ENSISID attribute in the Group by list.

iii) Click on Start Processing.

After each inventory has been loaded and intersected against each SSSI polygon a table will be produced. The table contains the intersect areas, a separate column for each inventory and a separate row for each ENSIS ID.

Technical

The has been written in C++ using Qt Creator IDE and compiled using the Visual Studio 2010 Qt libraries.

The interface is written using the Qt libraries, a cross-platform library meaning that software can be re-compiled to run on Linux and Mac OS. However, the raster processing is currently undertaken using direct calls to the Windows API, which means that it cannot be compiled directly to run on other operating systems. The software will however run directly on Linux using Wine.

The software was written using the non-commercial version of Qt and is made available under the GNU General Public Licence (GPL) v3.

ShapeLib

The software uses Shapefile C Library V1.2 to read from shapefiles. ShapeLib is available under the GNU Library General Public Licence, which has been superseded by the GNU Lesser GPL.

MITAB library

The software uses the MITAB library for reading from MapInfo tables and MID/MIF files. The MITAB library is released under a bespoke open source licence.

How the software works

The software works by creating rasterised bitmaps of the polygon layers, intersecting these and then counting the resulting bits. It uses the Window GDI calls to perform these operations, which results in much faster performance than if converting the polygons to byte arrays was coded or vector processing was used.

// Create device context, bitmap and brushes

m_hdc = CreateCompatibleDC(NULL);

m_hbrush = CreateSolidBrush(RGB(255,255,255));

m_hpen = CreatePen(PS_NULL, 0, RGB(255,255,255));

m_hbitmap = CreateDIBSection(m_hdc, &bmi, DIB_RGB_COLORS, &m_pData, NULL, 0);

// Draw polygon

Polygon(m_hdc, aPoints, nVertices); // Use PolyPolygon for multipart polygons

// Intersect two bitmaps

BitBlt(m_hdc,0,0,m_nWidth, m_nHeight, pRasterBitmap->m_hdc, 0,0, SRCAND);

// Extract bits

BYTE byte = ((BYTE)m_pData + i + j * bitmap.bmWidthBytes);

The original version of the software created two large device-dependent bitmaps (700km x 700km) when the program was first run. These used kernel memory (which is a limited resource) and the largest size possible was 25m resolution raster.

The current version uses tiled bitmaps, 5000 x 5000 bits in size, which are created as required (i.e. when drawing a polygon to that part of the overall image). It uses device independent bitmaps (DIB), whose memory can be accessed directly. Where a polygon covers multiple bitmaps it must be drawn multiple times, and the vertices shifted accordingly.

The pen used to draw the outer edge of the polygon, uses a null style (PS_NULL). The produces accurate results when compared to a white (which underestimates) or a black pen (which overestimates). The null pen appears to colour in a bit if the polygon covers more than 50% of the bit, i.e. errors are averaged out.

Accuracy

The following example shows the variation in accuracy with raster resolution. Nb. In some cases the lower resolution provides closer results to vector processing than that of higher resolution processing due to the effect of errors being cancelled out.

Projections

The software assumes that the spatial files are both using the same projection and that coordinates are in metres. The software no longer has any restrictions on the overall extent of poylgons (subject to available memory), and can for example easily handle polygons covering the whole of England using a 10m raster.

Performance

The software currently performs twenty times faster than FME for equivalent intersect queries. The area figures results are typically within 1% of those produced by FME. In fact the software didn’t replicate the spurious results coming from FME, that presumably originate from geometry errors in the input spatial data.

The software may return an error a) when it runs out of memory or b) when there is an error in the spatial file being imported. The software is only likely to run out of memory if the raster resolution is set below 10m or the tool is used on larger areas than England. Errors in the spatial files created during conversion, which result in spurious coordinate values, will result in the tool returning an error.

The tool has been optimised to improve performance:

• When grouping by polygons the overlap is only calculated when the group by attribute changes.

• The counting of bits and the BitBlt calls are restricted to the extent of the polygons that have been loaded.

• When counting of bits, the data is extracted in four byte words and the word cases (all bits or no bits set) are tested first.

• The tiles are now stored in an indexed container, which means the relevant tiles can be accessed more quickly.

• A spatial index of intersect polygons are stored on first pass so only those intersecting input polygons need to be loaded when processing subsequent input layers.

Source code

The key files are:

Version history

Version 0.5

• Added option to group by to ‘Sum’ intersect areas or to only include the ‘First’ intersect area found.

• Added layer count to progress status text, e.g. (1/22).

• Added warning when empty polygons are drawn, which may indicate a corrupt shapefile

• Extent of intersect layers are stored on first pass and extent of intersect polygons are indexed to 10km from second layer onwards. Only intersect polygons whose bounding box intersects that of an input polygon are loaded.

• Added separate (modal) progress bar with cancel button.

• Added option to intersect input polygons as well as to merge them.

Version 0.4

• Correctly releases memory for each SHPObject after retrieving row from shapefile.

• When calculating area, only retrieves bits within the overall extent of polygons that have been drawn.

• When intersecting bitmaps, limits the extent of the BitBlt operation to cover the polygons that have been drawn.

• Fixed crash when cancelling file open dialog.

• Simplified polygons after loading from spatial file (removing duplicate points when scaled to grid size).

• Display name of file being loaded and current row number during processing.

• When grouping, all adjacent polygons with the same attribute value are drawn before the area is calculated.

• Correct height and width are now used in call to BitBlt.

• Function to count number of bits set now retrieves data as a four byte word and first checks for ‘worst’ case, all bits set.

• Add call to GdiFlush prior to counting bits (GetExtent) to ensure that GPU in sync when processing very large number of polygons.

• Added column to results with the area of the polygon used to intersect the input polygon

• Clear button now removes individual spatial file name if selected, or clears all if none selected.

• Added ability for user to set the raster size, the default is 10m

• Added support for shapefiles containing polygonm and polygonz types.

• Tiled bitmaps are now stored in an indexed table (QMap) for faster lookup

• Tiled bitmaps are now created as required with no geographical limit on the overall extent.

• Properties form now shows overall extent of last processed input layer(s).

Version 0.3 (10th April 2011)

• Support for tiled bitmaps covering up to 700km x 700km at 10m resolution.

• Ability to calculate intersect areas against individual polygons in one of the layers and group results by one of the attributes.

• Ability to produce separate results for each of the input layers.

• Results output to a text file in tab separated table of areas for each layer versus group by attribute.

Version 0.2 (17th August 2011)

Support for two bitmaps 700km x 700km at up to 25m resolution. Ability to add multiple layers (shapefiles/MapInfo tables) to each layer (counting the combined area), intersect the polygons and calculate the intersect area.

Further development

At present multiple layers can be handled separately, but it is not possible to process a input spatial file containing multiple layers. This requires the ability to group by the input layer. This would need to optimised for repeated polygons of the same layer, otherwise it would probably be more efficient to do vector-based processing.

The software currently doesn’t support ESRI file-based geodatabases. ESRI have produced an API, but this only works with ArcGIS 10 geodatabases. It should be possible to read geodatabases using ArcObjects, although this would require ArcGIS and the development IDE (currently Qt Creator) to be installed on the same machine.