Extracting contours

One of my contour datasets has areas that are covered by a 10m contour interval whereas the majority is only covered by a 20m contour interval.

Contours - 10m versus 20m

So if I want to make a consistent 20m DEM I need to remove these pesky 10m contour intervals.

QGIS 1.9 > Select features using an expression

Select using expression

The expression used requires some regex (^[-+]?[0-9]*\\\.?[0-9]+$) to find any contours that have an altitude value that becomes a floating point number when divided by 100 and then halved.

The contours can then be deleted or the expression changed to regexp_match(((“ALTITUDE”/10)/2), ‘^[-+]?[0-9]*\\\.?[0-9]+$’) = 1 to select the desired contours and saved out as a new file.

Sources: http://www.regular-expressions.info/floatingpoint.html, http://gskinner.com/RegExr/

Advertisements

Automating geometry_columns extents

I think this sets a record for becoming irrelevant in the shortest space of time… the most recent QGIS 2.0 dev-build has been changed so that it automtically calculates the extents of datasets in the geometry_columns table without the need for the extra columns. The dangers of dealing with actively developed open source!

Having populated min_x, min_y, max_x and max_y columns (all int/null) in the geometry_columns table in SQL Server (automatically created when using ogr2ogr) creates huge time savings when loading SQL Server tables into QGIS 2.0.

To add these column to geometry_columns, use the following SQL:

ALTER geometry_columns
ADD min_x int null,
min_y int null,
max_x int null,
max_y int null

The stored procedure below automatically populates the extent values of all records in the geometry_columns table. It is currently set to only populate records with a particular SRID. Either change the SRID to the appropriate value or remove altogether:

CREATE PROCEDURE usp_extentsUpdate

AS

DECLARE @tablename VARCHAR(MAX)
DECLARE @srid INT = 28355

BEGIN

SET NOCOUNT ON;
DECLARE extents_cursor CURSOR FOR
SELECT f_table_name
FROM geometry_columns
WHERE srid = @srid

OPEN extents_cursor;

FETCH NEXT FROM extents_cursor
INTO @tablename

WHILE @@FETCH_STATUS = 0
BEGIN
DECLARE @sql VARCHAR(MAX)
SET @sql = ‘update geometry_columns set min_x = a.minx, min_y = a.miny, max_x = a.maxx, max_y = a.maxy from (select ‘ + ”” + @tablename + ”” + ‘ as table_name, MIN(ogr_geometry.STPointN(1).STX) as minx, MIN(ogr_geometry.STPointN(1).STY) as miny, MAX(ogr_geometry.STPointN(1).STX) as maxx, MAX(ogr_geometry.STPointN(1).STY) as maxy from ‘ + @tablename + ‘) as a where f_table_name = a.table_name’

EXEC (@sql)

SELECT @tablename + ‘ completed’

FETCH NEXT FROM extents_cursor
INTO @tablename
END
CLOSE extents_cursor
DEALLOCATE extents_cursor
END

QGIS – using contours to create a shaded relief

Creating a relief map is not something I have to do very often. Terrain rarely changes so unless the geographic scope or quality of the elevation datasets improve there is no real call to create relief data on a regular basis.

So when I needed to create a new shaded relief raster the other day it came as something of a shock that I couldn’t remember the process from the last time I had done it. Enter Google search engine and QGIS…

The easiest means of using QGIS to create a new shaded relief raster dataset from contours is through the GRASS (Geographic Resource Analysis Support System) plugin which is included in the QGIS install (using the Windows standalone installer) and the DEM relief shader plugin.

Step 1: create new GRASS mapset (Plugins > GRASS > New mapset)

The explanations about the data structure for a GRASS mapset are included in each step of the process but in short the ‘folder structure’ works like:

  • Database
    • Location
      • User 1 mapset
      • User 2 mapset
      • etc

Step 2 – add vector contour layer to QGIS (Layer > Add vector layer)

Step 3 – import vector layer into GRASS (Plugins > GRASS > Open GRASS tools > Modules List)

Use v.in.ogr.qgis to import the QGIS contour vector dataset into your GRASS mapset.

Step 4 – edit GRASS settings  (Plugins > GRASS > Edit Current GRASS Region)

It is important to ensure the geographical region you are working with is defined correctly so that any data processing is only confined to the area of interest and, therefore, speeding things up. Unless you happen to know the extent of your dataset and can type it in, click and drag on the QGIS canvas to define your area of interest.

Change the resolution so that the cell width and height match the contour interval of your data.

Step 5 – convert vector contours into raster (Plugins > GRASS > Open GRASS tools > Modules List)

Use v.to.rast.attr and set the attribute field to the elevation.

raster_contours

Raster contours

Step 6 – create a surface from the raster contours  (Plugins > GRASS > Open GRASS tools > Modules List)

Use r.surf.contour to create a surface. When the surface fist appears it will look pretty much like a large grey square. Select the layer then go to Layer > Properties > Style > Single band properties > Color map and change this to ‘Pseudocolor’ to get a better visual effect.

Surface

Surface

Step 7 – export surface to a geotiff  (Plugins > GRASS > Open GRASS tools > Modules List)

Use r.out.gdal.gtiff to output the surface to a geotiff that can be manipulated in QGIS. Choose the format based on the range of your data. I used UInt16 as this has a range of values from 0 to 65,535, there is no areas below sea level so a lower limit of 0 is fine and the data type Byte does not go high enough for some of the mountains in the area of interest.

Step 8 – create a shaded relief (Plugins > Shaded Relief > Shaded Relief)

Select the geotiff surface in the layers list to make it the active layer before opening the DEM relief shader tool. Simply leave all the settings as per the defualts and click OK. The output file will be written to the same location as the surface file.

Shaded relief

Shaded relief

And there is your relief map!

Useful References:

http://linfiniti.com/2010/12/3d-visualisation-and-dem-creation-in-qgis-with-the-grass-plugin/