Spatially index (*.qix) all shapefiles in a folder (and sub-folders)

pushd D:\<folder>

for /R %%f in (*.shp) do ogrinfo -sql “create spatial index on %%~nf” “%%f”

Point the function at the correct folder. Loop recursively through the folder/sub-folder structure looking for shapefiles then undertake an ogrinfo SQL command on each file. %%nf provides the file name and %%f provides the full path to the file. Double ‘%%’ is required when the commands are run from a batch file, otherwise single ‘%’ will do.

Splitting delimited text in SQL

Split delimited text:

declare @userinput varchar(100) = ‘GREAT ALPINE ROAD – BRIGHT’
declare @roadname varchar(100)
declare @locality varchar(100)

select
@roadname = rtrim(PARSENAME(replace(@userinput, ‘-‘,’.’),2)),
@locality = ltrim(PARSENAME(replace(@userinput, ‘-‘,’.’),1))

declare @geo geometry

select @geo = ogr_geometry
from [GISProduction].[dbo].[locality_polygon]
where locality = @locality

select pfi
from [GISProduction].[dbo].[tr_road] as a
where a.ogr_geometry.STIntersects(@geo) = 1
and a.ezi_rdname = @roadname

Selecting the first record from grouped records

I need to create a table of localities with the number of waste bins associated with each, however, the address table can have multiple localities for a single property number. Therefore, rather than skew the results by counting a single bin twice (or more) for each locality the property has an address in it is better to simply select the first locality in the list of address and assign that to the bin count.

For example:

select b.property_address_2 as locality,

sum (a.number) as bin_count

from rd_waste as a

join

(select distinct property_number,

property_address_2,

row_number() over (partition by property_number order by property_address_2) as record 

from rd_address) as b

on a.property_number = b.property_number

where a.description like ‘%Waste%’

and b.record = 1

group by b.property_address_2

order by b.property_address_2

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/

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

Find my nearest…

I currently use the following method for calculating the nearest ‘thing’ to a point of interest:  http://blogs.msdn.com/b/isaac/archive/2008/10/23/nearest-neighbors.aspx, my thanks to Isaac Kunen!

With some slight modifications…

This stored procedure calculates the distance from the centroid of a parcel to the closest fire plug, an important measurement when subdividing land to ensure all part of a new property can be reached by a fire hose!

Some notes:

  • @XPoint and @YPoint are parcel centroid parameters passed to the stored procedure in the form of UTM coordinates
  • When creating the geometry using the point parameters it is important to get the SRID correct, this caused a minor headache for a while until I figured out that I needed to use 28355 (Zone 55)
  • A full explanation of the SQL used here can be found in the referenced web page as well as the SQL needed to create the ‘numbers’ table that is used to iterate the search area
  • This SQL is used as part of our Online Property Report mapping service: http://maps.alpineshire.vic.gov.au/opr/

CREATE PROCEDURE [dbo].[usp_GetFirePlugDist]

@XPoint VARCHAR(100),
@YPoint VARCHAR(100)

AS

DECLARE @Point VARCHAR(100)
DECLARE @start FLOAT = 1000;
DECLARE @x GEOMETRY;

SET @Point = ‘POINT(‘ + @Xpoint + ‘ ‘ + @Ypoint + ‘)’

SET @x = GEOMETRY::STGeomFromText(@Point, 28355);

WITH NearestPoints AS
(
SELECT TOP(1) WITH TIES *, Round(b.ogr_geometry.STDistance(@x), 0) AS DISTANCE
FROM numbers AS a
JOIN fire_plugs AS b
ON b.ogr_geometry.STDistance(@x) < @start*POWER(2,a.n)
ORDER BY a.n
)
SELECT DISTINCT TOP(1) DISTANCE AS distance_metres
FROM NearestPoints

GO

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/