Postgres Command Line : psql : Create a switch to determine whether an Arc has a clockwise or anti-clockwise direction

This is a post focused around Network distance measurement using POSTGIS and routing plugins.. This discusses the use of Ordnance survey Road Network information and its preparation in using it for routing.

The Ordnance Survey open source road network layer includes roundabouts that have an attribute value of roundabout. Great but looking at them closely some of the constituent arcs are digitised in a clockwise direction while others are digitised in an anti-clockwise direction. When using dijkstra routing with weighting to prevent incorrect pathing it is necessary to ensure that networks are weighted in the correct fashion. Directional weighting only works if you know the direction and you weight appropriately. For use with directional routing in the UK roundabouts any directional weighting should prevent travel in anticlockwise direction. ST_reverse will correct incorrect direction BUT the Ordnance survey layer seems to have no attribute that consistently indicates whether an arc on a roundabout has or has not been digitised correctly. Marking lines with direction and viewing them we see the random nature of many arcs on roundabouts.

Here is Sheriff Hall Roundabout on Edinburgh City Bypass

Here is Straiton roundabout just north of the bypass

and finally Hermiston Gate roundabout again following the theme on Edinburgh city bypass

It got me thinking was there a way to determine whether arcs on roundabouts could be determined to be clockwise or anti-clockwise?

Having thought about it in my head quite a bit I determined that it probably would be possible if we knew three points and could create some kind of virtual graph with the start point being at 6 and a finish at 12 with mid points at 9 if travelling in a clockwise position and 3 if travelling in an anti-clockwise position.

I had a look around and the following post seemed to tally with the idea of three points and positive and negative relating to clockwise or anticlockwise.

Maths to determine direction of an arc

Having looked at this I set about working through the problem in Excel to see if I could get consistent results.

Firstly I created a set of directions North West South and East and placed some coordinates that could be used in calculations.

I then went forward and tested whether I could identify the direction of various arcs from these coordinates in excel using the formula identified on Stack Exchange.

Here I replaced a,b,c with b Beginning, m Middle and f Finish

And I decided to work through manually in excel to ensure that I had the linear algebra correct.

and further testing figures

So firstly I create a separate table that just shows the roundabouts

CREATE TABLE t001roundaboutsntroadlink as select id pkidt001, st_linemerge(geom) as geom from ntroadlink where formofway = 'Roundabout';

In the above I use st_linemerge to ensure that all geometry is linestring as this is necessary to use the st_startpoint and st_endpoint postgis functions.

Next I added the the required variables from stack overflow algebra to the line table

ALTER TABLE t001roundaboutsntroadlink add column bx float(8),
Add column by float(8),
Add column mx float(8),
Add column my float(8),
Add column fx float(8),
Add column fy float(8),
Add column ux float(8),
Add column uy float(8),
Add column vx float(8),
Add column vy float(8),
Add column uxtimesvy float(8),
Add column uytimesvx float(8),
Add column uxv float(8);

Next I needed to identify a b beginning, m middle and f finish point for each line that I wanted to test.

b points (beginning)

CREATE TABLE t002bpoints AS SELECT pkidt001 as pkidt001,st_startpoint(geom) as geom, st_x(st_startpoint(geom)) as bx, st_y(st_startpoint(geom)) as by from t001roundaboutsntroadlink;

m points (middle)

CREATE TABLE t002mpoints AS SELECT pkidt001 as pkidt001,st_lineinterpolatepoint(geom,0.5) as geom, st_x(st_lineinterpolatepoint(geom,0.5)) as mx, st_y(st_lineinterpolatepoint(geom,0.5)) as my from t001roundaboutsntroadlink;

f points (finish)

CREATE TABLE t002fpoints AS SELECT pkidt001 as pkidt001,st_endpoint(geom) as geom, st_x(st_endpoint(geom)) as fx, st_y(st_endpoint(geom)) as fy from t001roundaboutsntroadlink;

It was then a case of simple update queries to complete the table

update t001roundaboutsntroadlink set bx = st_x(st_startpoint(geom));
update t001roundaboutsntroadlink set by = st_y(st_startpoint(geom));
update t001roundaboutsntroadlink set mx = st_x(st_lineinterpolatepoint(geom,0.5));
update t001roundaboutsntroadlink set my = st_y(st_lineinterpolatepoint(geom,0.5));
update t001roundaboutsntroadlink set fx = st_x(st_endpoint(geom));
update t001roundaboutsntroadlink set fy = st_y(st_endpoint(geom));
update t001roundaboutsntroadlink set ux=mx-bx;
update t001roundaboutsntroadlink set uy=my-by;
update t001roundaboutsntroadlink set vx=fx-mx;
update t001roundaboutsntroadlink set vy=fy-my;
update t001roundaboutsntroadlink set uxtimesvy = ux*vy;
update t001roundaboutsntroadlink set uytimesvx= uy*vx;
update t001roundaboutsntroadlink set uxv = uxtimesvy-uytimesvx;

Labelling up the roundabouts Hermiston Gate now looks like

And Sheriff Hall Roundabout now looks like this

Compared with a correctly directed roundabout

CREATE TABLE t001roundaboutsntroadlinkcorrected AS TABLE t001roundaboutsntroadlink;

And now correct the items display as previous and see what we see.

UPDATE t001roundaboutsntroadlinkcorrected set geom = st_reverse(geom) where uxv > 0;

Sheriff hall roundabout now

and some proof that reasonable number of lines were updated.

Which is an indication that all roundabouts arcs have been corrected properly

But a zero uxv value indicates a straight line.

It should however be possible to match starts with finishes for overlying points and where a line has 0 value of uxv and its ends and finishes are not matched with adjacent opposites create a switch to reverse the direction of all lines that are incorrect compared to their neighbours thus only correcting incorrect directions. Haven’t done that in this case.

QGIS – Working with Free BASEMAP Services for Desktop Projects (Bing and Google Maps Aerial and Lines)

As soon as you need to practically implement any information in Spatial Databases display of the information through a mapping front end becomes absolutely vital. Many database administrators are used to simply going into their favourite database editors and displaying the raw subsets of tables and queries. That works well for financial transactions and inventory tables were collapsing the attributes of objects into single digits is often valuable or possibly preferable to simple photos. When dealing with boundary information the complete opposite applies. Display of information as simple screens of matrix numbers is completely useless.

Additionally often boundaries make zero sense unless referenced to the land beneath them either through aerial photography or topographic maps.

In a previous time where I worked we actually commissioned a company to give us aerial photos of a local authority. This was not an insignificant amount of money and was probably only marginally reduced by the vendor having multiple channels of sale. Google and Microsoft are now very good in offering very good aerial and line interpretations for limited use to companies and individuals. This is absolutely great as it can be used as background either to confirm accuracy of other information or as data upon which to calculate further information (eg routing).

So how can an individual get up and started with some of these basemaps.

Sometime recently (I know not when exactly) QGIS changed its implementation of Open Street Maps through their desktop – rather than being an additional plugin Open Street Map provision is now included on install.

Here I am working with QGIS version 3.10

Now you should be presented with the Data Source Manager Dialog which looks like this

Now expand the XYZ Tiles

Open a project in QGIS and then select OpenStreetMap and select Add Layer to Project

And mapping should appear.

And of course Open Street Map has global coverage so this should work for anyone anywhere in the world.

So there are other sources of basemap information

The following dialog should then appear.

And here are important urls supplied by Google / Microsoft that are very useful.

Google Hybrid
{x}&y={y}&z={z}

Google Satellite
{x}&y={y}&z={z}

Google Road / Streets
{x}&y={y}&z={z}

Bing Aerial
{q}.jpeg?g=1

Thank you to Google and Microsoft

021 Postgres with PostGIS plugin – Create junction table sites in catchments

Quick post that I will come back and edit

So we need two tables

t001asites which has a geometry field called geom
and another table which will be the catchments table called
t002bcatchments which has a geometry field called geom.

Both tables must have a serial primary key of pkid and both tables must be polygon data and the geom field MUST be defined as polygon and NOT multipolygon.

Air code is as follows.

    1. Create table containing digitised polygons of housing sites.
    2. Create table containing digitised polygons of catchments.
    3. Measure the area of the housing sites and place that value in an area column within the housing sites table t001asites.
    4. Split the housing sites by the catchment boundaries ensuing that each split polygon inherits the catchment it was split by.
    5. Re-measure the areas of these split sites and add an area column to store the new calculations.
    6. Divide figure obtained in 5. by figure obtained in 3 which will indicate the proportion of the housing site is in which catchment.
    7. Perform a least remainder method on the individual sites grouped by their original housing sites to ensure the proportions sum to 1.

So to the code

BEGIN;
SET LOCAL check_function_bodies TO FALSE;
CREATE OR REPLACE FUNCTION part01catchjunctionmaker() returns void as $$
Alter table t001asites add column area integer;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part02catchjunctionmaker() returns void as $$
Update t001asites set area=ST_Area(geom);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part022catchjunctionmaker() RETURNS void AS $$
DROP TABLE IF EXISTS t200;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part03catchjunctionmaker() RETURNS void AS $$
CREATE TABLE t200 AS select a.pkid as t001pkid, b.pkid as t002pkid, a.area as t001area, ST_intersection(a.geom, b.geom) as geom FROM t001asites a, t002bcatchments b;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part04catchjunctionmaker() RETURNS void AS $$
ALTER TABLE t200 add column pkid serial primary key, add column area integer,add column proportion decimal (10,9);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part06catchjunctionmaker() RETURNS void AS $$
UPDATE t200 SET area=ST_Area(geom);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part07catchjunctionmaker() RETURNS void AS $$
DELETE from t200 where area=0 or null;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part08catchjunctionmaker() RETURNS void AS $$
UPDATE t200 SET proportion= cast(area as decimal)/cast(t001area as decimal) WHERE area > 0;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part088catchjunctionmaker() RETURNS void AS $$
DROP table IF EXISTS t201;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part09catchjunctionmaker() RETURNS void AS $$
Create table t201 as Select pkid,t001pkid,t002pkid, t001area, area, proportion, sum(proportion) OVER (PARTITION BY t001pkid ORDER BY t001pkid, proportion) as cum_proportion FROM t200 ORDER BY t001pkid, proportion;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part10catchjunctionmaker() RETURNS void AS $$
Alter table t201 add column value decimal (14,9),
Add column valuerounded integer,
Add column cumulvaluerounded integer,
Add column prevbaseline integer,
Add column roundproportion integer;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part11catchjunctionmaker() RETURNS void AS $$
UPDATE t201 set value = proportion * 100;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part12catchjunctionmaker() RETURNS void AS $$
UPDATE t201 set valuerounded = round(value,0);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part13catchjunctionmaker() RETURNS void AS $$
update t201 set cumulvaluerounded = round((cum_proportion*100),0);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part14catchjunctionmaker() RETURNS void AS $$
update t201 set cumulvaluerounded=100 where cumulvaluerounded = 101;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part15catchjunctionmaker() RETURNS void AS $$
update t201 set prevbaseline = round((cum_proportion - proportion)*100);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part16catchjunctionmaker() RETURNS void AS $$
update t201 set roundproportion = (cumulvaluerounded-prevbaseline);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part17catchjunctionmaker() RETURNS void AS $$
DELETE from t201 where roundproportion=0 or null;
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part18catchjunctionmaker() RETURNS void AS $$
alter table t201 add column proppercent decimal(3,2);
$$ LANGUAGE SQL;
CREATE OR REPLACE FUNCTION part19catchjunctionmaker() RETURNS void AS $$
update t201 set proppercent = cast(roundproportion as decimal)/100;
$$ LANGUAGE SQL;
COMMIT;

and now a function to pull it all together;

CREATE OR REPLACE FUNCTION createcjt()
RETURNS TEXT AS
$BODY$
BEGIN
PERFORM part01catchjunctionmaker();
PERFORM part02catchjunctionmaker();
PERFORM part022catchjunctionmaker();
PERFORM part03catchjunctionmaker();
PERFORM part04catchjunctionmaker();
PERFORM part06catchjunctionmaker();
PERFORM part07catchjunctionmaker();
PERFORM part08catchjunctionmaker();
PERFORM part088catchjunctionmaker();
PERFORM part09catchjunctionmaker();
PERFORM part10catchjunctionmaker();
PERFORM part11catchjunctionmaker();
PERFORM part12catchjunctionmaker();
PERFORM part13catchjunctionmaker();
PERFORM part14catchjunctionmaker();
PERFORM part15catchjunctionmaker();
PERFORM part16catchjunctionmaker();
PERFORM part17catchjunctionmaker();
PERFORM part18catchjunctionmaker();
PERFORM part19catchjunctionmaker();
RETURN 'process end';
END;
$BODY$
LANGUAGE plpgsql;

016 Postgres command line : psql : Strip out the Z coordinate from a geometry field

When creating a topology the geometry field cannot contain a Z coordinate.

OK but the Ordnance Survey Open Data highways layers containse a Z coordinate. Previously I had stripped this out using the latest version of QGIS which has a tick box in the front end that allows for import stripping of the z coordinate in the process. If you don’t have access to the latest QGIS version how can you strip out the z coordinates.

ST_FORCE2D

ALTER TABLE public.nuroadlink ADD COLUMN geom2(multilinestring,27700);
UPDATE public.nuroadlink SET geom2 = ST_FORCE2D(public.nuroadlink.geom);
ALTER TABLE public.nuroadlink drop column geom;
ALTER TABLE public.nuroadlink RENAME COLUMN geom2 TO geom;

013 Postgres command line : psql : Using ST_Within function to build junction tables to compare 2 separate polygon tables

First off let us create a new database to hold our examples in.

CREATE DATABASE stwithindb;

Now add the postgis extension.

Lets create two tables one called fields and one called plots

CREATE TABLE
t00001fields
(
pkid serial primary key,
fieldname varchar(50),
geom geometry(polygon,27700)
)
;
CREATE TABLE
t00002Plots
(
pkid serial primary key,
plotname varchar(50),
geom geometry(polygon,27700)
)
;

Now lets go to QGIS connect to the PostGIS instance add the tables and create some test data manually.

Here I have added fields in green with bold number labels and plots in brown with smaller number labelling. The numbers represent the pkid fields.

Now here I can quickly run a query to identify the plots that are in fields

SELECT t00002plots.pkid
FROM
t00002plots,
t00001fields
WHERE 
ST_WITHIN(PUBLIC.T00002PLOTS.GEOM, PUBLIC.T00001FIELDS.GEOM);

And it correctly identifies that plot 1 is within the fields layer.

But what would be great in an application is to have some kind of junction table that individual master records could display their children on. For this we need a junction table that links between the field and plots table showing the pkids from each.

SELECT t00002plots.pkid as Plotspkid,t00001fields.pkid as Fieldspkid
FROM
t00002plots,
t00001fields
WHERE 
ST_WITHIN(PUBLIC.T00002PLOTS.GEOM, PUBLIC.T00001FIELDS.GEOM);

Now I will move plot 2 into field 3 and rerun the above.

The layer now looks like

and running the former query we get.

Now its possible to either create a junction table to hold this information..

eg

CREATE TABLE t00010fieldplotjunction AS 
SELECT t00002plots.pkid as Plotspkid,t00001fields.pkid as Fieldspkid
FROM
t00002plots,
t00001fields
WHERE 
ST_WITHIN(PUBLIC.T00002PLOTS.GEOM, PUBLIC.T00001FIELDS.GEOM);

or we can create a view that will constantly calculate this everytime it is seen

CREATE VIEW v001FieldPlotJunction AS
SELECT t00002plots.pkid as Plotspkid,t00001fields.pkid as Fieldspkid
FROM
t00002plots,
t00001fields
WHERE 
ST_WITHIN(PUBLIC.T00002PLOTS.GEOM, PUBLIC.T00001FIELDS.GEOM);

Now if I add a few more plots and fields and then pull up the view we shall see that everything has been adjusted

and running the view we now get

In some circumstances this calculation may be expensive so we may wish to run and create a junction table overnight other times we may be happy to do it fully dynamically. Of course in a front end you could query and filter such that only one record was compared against the fields plot at anytime. Very useful nonetheless.

011 : Postgres amalgamate consecutive lines into a single line in a table

Here we take much of the work covered in post 010 and take the parts and user st_union to merge into a single record and place it in a table created by transforming a view into a table

Firstly go to your psql line and ensure that you are logged in with a username that you wish to be the owner of the table. In my case general

logging into edinburgh routing database

Now same measurement as before but this time we shall make a view out of the measurements then load that into a new table before deleting the view leaving us with the table with a combined measurement.

CREATE VIEW v001firstmeasurement AS SELECT seq,  id1 AS node, id2 AS edge, cost, geom, agg
  FROM pgr_dijkstra( 'SELECT id, source, target, st_length(geom) as cost FROM public.t01roadnetwork', 15883, 10967, false, false  ) as di
  JOIN public.t01roadnetwork pt ON di.id2 = pt.id ;

CREATE TABLE t003 as select sum(cost), st_union(geom) from v001firstmeasurement;

DROP VIEW v001firstmeasurement;

It is important in notepad to remove the blank spaces in the editor this looks as follows.

We then should then get some kind of confirmation that the view and table are created before the view is then dropped again. There might be a more efficient way of doing this but this was my first experiment. And we can go back to QGIS 3.4 and display the now single line in our project. Complete with now accurate measurement. It should be noted that if you were wanting to do multiple line measurements you would need to step out of the create statement and use an insert statement for all subsequent insertions as follows.
insert into t003(sum,st_union) select sum(cost),st_union(geom) from v001firstmeasurement;
This would allow you to do multiple measurments. I haven’t added up the measurement but it looks about right.

010 Postgres command line : psql : Getting started with pgrouting using open data from Ordnance Survey to identify and measure the shortest route between two points.

Objective here is to write a series of queries that can be used to measure the shortest distance between selected paired locations on a network such that the geometry of the routes can be calculated and displayed on a map.

For this particular tutorial you will need – QGIS 3 or higher and a version of Postgres I am using version 11.0 here (I have upgraded since my former posts). I believe this tutorial will work with previous versions but if you are following along now might be a good time to upgrade.

QGIS 3.4 or higher – needed as the Ordnance Survey road network geometry contains a z coordinate which will prevent the creation of the required geometry for measurement. QGIS 3 introduced the ability to save geometry excluding z coordinate. If you have a network without z coordinates you should not require this.

So let us first get the data. Here you tick the option in the top right hand corner – scroll to the bottom and submit your request after which you will be asked a few basic questions along with email address you wish the download to be sent to after a few minutes you should be sent the download link through your email – follow the instructions and you should be able to get the information

Ordnance Survey Open Data

The information you are downloading is a block framework for the whole of the uk. When you unzip the download into a folder you will see multiple files. We will be using a section of the national dataset relating to Edinburgh – NT. Choose the block or selection that you are interested in. More blocks may take more time however.

Open QGIS
Create a new project : eg EdinburghRouting.qgz
Load in your chosen network block : eg NT_RoadLink.shp

Select the layer you just loaded in : eg NT_RoadLink.shp

and navigate to the following in the menu settings
Layer / Save As

Fill out the Save Vector Layer as … dialog box
IMPORTANT – ensure within the Geometry section
Geometry type is set to LineString
Include z-dimension is unticked

Give the new file a name : eg ntosroutingnetwork.shp

Hit ok

Within the layer dialog of QGIS your new layer should appear you can now remove the for NT_RoadLink shape file from the project

Next go to your version of PostgreSQL and using a superuser account create a new database : eg edinburghrouting

I would suggest you use lower casing as well

As a superuser ensure you add the postgis and pgrouting extensions.

Next I set up the following connection between the QGIS project and PostgreSQL

Personal tastes may vary but I like like to select
Also list tables with no geometry
Allow saving/loading QGIS projects in the database

OK the selection and you should now have a connection to the new database you just created.

QGIS has an excellent dbmanager window which we will use to load our new shape file which excludes the z layer into the new database we created in PostgreSQL

Ensuring that you have a connection to your localpostgis database hit the

ImportLayerFile

Here I load the information into a new table t01roadnetwork

On pressing OK there will be delay after which if things go well you will receive the following message.

As ever it is good to check that things appear to be going well.
Add the layer to your project and determine approximately whether import was successful.

Next back in psql command line and in an editor we are going to run 4 queries
The first 2 add columns that are required in the shortest distance algorithm we shall use, the third will allow anyone to write an aggregation function to see the total cost of the route and the last creates a topology for the road network.

alter table public.t01roadnetwork add column source integer;
alter table public.t01roadnetwork add column target integer;
alter table public.t01roadnetwork add column agg smallint default 1;
select pgr_createTopology('public.t01roadnetwork', 0.0001, 'geom', 'id');

If things go correctly you should see the database engine start to create the topology and what I see is it gradually stepping through the creation process.

and on completion you should have something like the following:

A new table has been added to the edinburghrouting database and next step is to display the network and its vertices. In QGIS.

In QGIS we should see something like

The next thing that I like to do is to label the nodes so that for quick identification.

And look to the t01roadnetwork table and see if the columns are clear and present.

We are now ready to make a measurement. Here I choose the nodes 15883 and 10967

SELECT seq, id1 AS node, id2 AS edge, cost, geom , agg
  FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(geom) as cost FROM public.t01roadnetwork',
    15883, 10967, false, false
  ) as di
  JOIN public.t01roadnetwork pt
  ON di.id2 = pt.id ;

Now we can load this as a new layer and then improve the symbology

Doing this we get.

It should be noted that the line you see is a collection of lines. In my next post I will go through and indicate how we can amalgamate that into a single line for storage in a table.

Congratulations if you have got this far you should be able to measure the shortest distance between any two points on a valid network by altering the numbers.

Add Open Street Map to Background QGIS Project and then Digitise against imported Raster

The following is a workflow that can be used to get a raster basemap of anything into QGIS which you then reference to Open Street Map Layers ready for digitising against. This will be useful for approximate digitising of masterplans and approsimate digitisation of housing completions.

Firstly ensure you have dowloaded QGIS and added the following two plugins
OpenLayers Plugin

Georeferencer GDAL
Plugin

Opening QGIS now lets add the the Open Street Map Raster

From a blank project selection of Open Street Map should give you the following result

Now zoom to the approximate location where you wish to have a unique basemap. You will be referencing points on this map to points on your imported raster so you should zoom into a location to the extent that you can identify common locations between the two maps.

Identify the basemap you wish to have in your particular QGIS map here I choose freely available masterplan from Calderwood development in West Lothian from planning application 0524/P/09

Within the menus navigate to
Raster / Georeferencer /

You should be presented with the following window.

Hit the add raster button in the top left

Select the basemap you wish to add to your project and ensure that the coordinate system that you choose is OSGB 1936 / British National Grid

Next you want to add reference points to the basemap that will allow for you to put the basemap against it – This is done using the button marked

Next hit the settings button

You should now be presented with the Transformation parameters windows dialog as follows.
The dialog will remember old parameters if not ensure that you have the same selections (with your own selection of output raster location) as mine.

Now hit the play button the raster will be added to your map and the georeferencer will be reduced and moved to the bottom left of the corner where you will be open it and reduce it in size if you wish. You can now go in and alter the transparency so that it is possible to see both Open Street Map and your newly added raster

You should now be presented with something like the following – if there are red dots on the screen this is because you have not closed georeferencer down – simply open the window up again and hit file close.

VBA Function to Create Table of Import strings using OGR2OGR targeting a SQL Server

Do you have many shape files you wish to import into a local SQL Server Database so that you can display them in QGIS or serve them on Geoserver?
Here’s a short function I wrote that will take a table called T0001OpenStreetMapLayers with fields PKID/Name/Directory/Type/Flag – and produce OGR strings that can then be used to load them into a local SQL Server / SQL Express or SQL Azure

For this to be useful you will need
A version of QGIS
A local SQL Server copy (in this case SQL Server Express)
A database within your copy called OpenStreetMap
All shape files in the same directory
You will also need to figure out how to get all those shape files into the table T0001OpenStreetMapLayers table
A starting database with 2 tables
T0001OpenStreetMapLayers with populated fields PKID/Name/Directory/Type/Flag
T0002OGRStrings blank table with fields PKID/CommandLine – This is where all the Command Line Strings will be stored

Public Function CreateTableOGR2OGRString()

Dim rs1 As DAO.Recordset
Dim rs2 As DAO.Recordset
Dim db As DAO.Database
Dim O2O As String
Dim LCounter As Integer
Dim strQuote As String
Set db = CurrentDb
strQuote = Chr$(34)


LCounter = 1
While LCounter < 3000
LCounter = LCounter < 3000

Set rs1 = CurrentDb.OpenRecordset("SELECT T0001OpenStreetMapLayers.PKID, T0001OpenStreetMapLayers.Name, T0001OpenStreetMapLayers.Directory, T0001OpenStreetMapLayers.Type, T0001OpenStreetMapLayers.Flag FROM T0001OpenStreetMapLayers WHERE (((T0001OpenStreetMapLayers.Type)=1) AND ((T0001OpenStreetMapLayers.Flag)=0 Or (T0001OpenStreetMapLayers.Flag) Is Null));")
O2O = "ogr2ogr -append -f MSSQLSpatial " & strQuote & "MSSQL:server=DESKTOP-JECT7QO\SQLEXPRESS;database=OpenStreetMap;trusted_connection=yes" & strQuote & " " & strQuote & rs1!Directory & rs1!Name & ".shp" & strQuote & ""


rs1.Edit
rs1!Flag = 1
rs1.Update
rs1.MoveNext
rs1.Close

Set rs2 = CurrentDb.OpenRecordset("T0002OGRStrings")
With rs2
.AddNew
rs2!CommandLine = O2O
rs2.Update
rs2.Close
End With
Wend
End Function

For SQL Azure target databases replace the yellow connection string with something resembling;

MSSQL:Server=tcp:azureinstance1.database.windows.net;Database=TouristDB1;
Uid=tom@azureinstance1.database.windows.net;Pwd=Edinburgh;

There are multiple methods of finding the name of your SQL Instance – Ignoring the fact that you won’t be able to connect to it if you don’t know it – Within SSMS you can right click on the instance and look to properties but the name itself is usually in the instance path of SSMS as well.

Creation of SITE History from Planning Application Polygons using QGIS

In planning it is important to know the planning history on a site. The status and likelihood of approved permission will often relate to previous permissions. Many council planning systems do not specifically relate planning applications to each other and there may be situations where you would like to create such links. This is essentially an excercise in using spatial analysis to create the junction table to hold what are many to many relationships.

If your datasets are in any way large you will need to set aside a computer so that it can perform the calculations. When I first tried this the process took a weekend with queries running overnight.

Start by obtaining as many years of planning application polygons as you can. Here I use polygon files in shape format.

The polygon file or shape file should be in one file so if you need to merge the shape files you have together. I did this and the file I ended up with was

AllPlanningApplications.shp

Next – Delete all attribute fields EXCEPT the planning application number.

Next – Create a centroids file from AllPlanningApplications.shp I called mine
AllPlanningApplicationsCentroids.shp

The next series of iterations are about getting a unique set of polygons with which we can go forward and generate a set of SITEPKIDS that can be attached to the child records.

Step – Using AllPlanningApplications.shp ADD an additional field called area and populate it using QGIS $area calculation – save this file.

Step – this is where it becomes interesting – in most authorities there are a vast number of planning application boundaries that overlap. Performing a dissolve at this point would result in a large congealed set of polygons that could not clearly identify unique sites. Thus buffering the polygons down we can start to identify unique sites. This is particularly important where boundaries are completely contiguous to each other.

sites the buffering command is used within the geometry tools to try to separate adjacent overlapping and contiguous polygons.

Step ‐ Create two files from the AllPlanningApplications.shp one for polygons less than 4500 metres squared and one for more than or equal to 4500 metres squared. This is to allow for two differing buffering processing to be performed on each.

AllSmallLessthan4500PlanningApplications.shp

AllLargeGreaterthanequal4500PlanningApplications.shp

Now the 4500 is an empirical figure that was subjectively chosen there may be a better figure feel free to investigate.

The following 2 steps also introduce empirical figures for the buffering that can be altered as appropriate.

Step ‐ Take the file AllSmallLessthan4500PlanningApplications.shp and create a buffer polygon file of this with

boundaries of less than 2m lets call it

AllSmallLessthan4500PlanningApplicationsBufferMinus2.shp

Step ‐ Take the file AllLargeGreatethanequal4500PlanningApplications.shp and create a buffer polygon file with

boundaries of less than 20m lets call it

AllLargeGreaterthanequal4500PlanningApplicationsMinus20.shp

THIS NEXT STEP WILL TAKE SEVERAL HOURS IT MAY BE BEST TO DO EACH ONE OVERNIGHT

Step ‐ Perform dissolves on both of these new files ensuring that dissolve all is used names could be something like

Vector / Geoprocessing Tools / Dissolve /

Set input layer alternatively to the two above files and set Dissolve field to dissolve all.

Suggested file Names are

MultipartDissolvedPolygonLessthan4500PlanningApplicationsBufferMinus2.shp

MultipartDissolvedPolygonAllLargeGreaterthanequal4500PlanningApplicationsMinus20.shp

Step You should now have two shape files of a large multipart polygon you want to perform the multipart to single part operation now

Vector / Geometry Tools / Multipart to Single Part

Processing involved with this is typically quick and suggested names for these files are

DistinctPolygonsAllSmallLessthan4500PlanningApplicationsMinus2.shp

DistinctPolygonsAllLargeGreatethanEqual4500PlanningApplicationsMinus20.shp

Add area column and identify the largest polygon on the small files

Add area column and identify the smallest polygon are on the large files you may want to remember this.

Step ‐ perform merge on these two files to get

Vector / Data Management Tools / Merge

CombinedSmallandLargeDistinctPolygonsPlanningApplicationswithbuffering.shp

ONGOING investigation ‐ would Difference be better than dissolve on this and should the above files be put together before

Step ‐ perform dissolve

Vector / GeoprocessingTools / Dissolve

ensure that ‐‐Dissolve all‐‐ is selected

DissolvedCombinedSmallandLargeDistinctPolygonsPlanningApplicationswithbuffering.shp

Step now you want to split mutlipart to single

DistinctPolygonsAllPlanningApplications.shp

Step add field called SitePKID and populate it using $rownum command.

Step

Vector / Data Management Tools / Join Attributes by Location

Set Target Vector Layer as

AllPlanningApplicationsCentroids.shp

Set Join Vector Layer as

DistinctPolygonsAllPlanningApplications.shp

Ensure that Keep all records (including non‐matching target records are kept)

Output Shapefile suggestions

AllPlanningApplicationsCentroidswithSitePKID.shp

If there are centroids without Site PKIDs put them to the end and give them consecutive unique row numbers. The attribute file associated with AllPlanningApplicationsCentroidswithSitePKID.shp should now be a child table of the shape file DistinctPolygonsAllPlanningApplications.shp perform checks here to see if all centroids within a polygon defined by the distinct polygons have the same SitePKID and that it is matched by the SitePKID of the Parent shape file.

You should be able to do a join on the this file to get the PKID back into the very original file.

AllPlanningApplications.shp

Finally perform a dissolve on the corrected AllPlanningApplications.shp file but this time dissolve on the field

SitePKID

You can call this

DistinctCorrectedPolygonsAllPlanningApplications.shp

QED!!!!

QGIS – Import shape file into PostGIS Table

The following uses
QGIS 2.14.2 Essen and
PostGres 9.5

A number of local authorities have released information through the UK’s data government site. The following example uses a shape file obtained from Lichfield District Council – At 2nd of October 2016 this was available for download from the following link

Lichfield Planning Applications

Open up QGIS and add Lichfield’s planning application shape file
qgisessen2142

Now scan along the top menu and go to Database

Select the sub menu DB Manager and then DB Manager

dbmanager

The following windows dialog should appear

dbmanagerdialog

Expand the area on the left named PostGIS – any PostGIS instances that you have created should be visisble here. Note you will have to have the PostGIS server running. Then highlight the actual instance that would like to import information into.

In this case I use the instance LocalPostGres

dbmanagerdialog

Choose the third icon from the left.
dbmanagerimportlayerfile

It should be noted that the window on the right may or may not show the correct connection to the database on the right.

importdialog

Name the table you wish to create and then hit OK – additional parameters are available.
There will be a delay before a confirmation of successful import happens – try to not issue commands during this time – once confirmation has been received go back into the PostGIS option and add the layer.

QGIS 2.8.1 – Useful Functions and Operators – Field Calculator

Calculate eastings and northings of centroid within polygon layer
xmin(centroid($geometry))
ymin(centroid($geometry))

Calculate area and perimeter of a polygon layer
$area
$perimeter

Calculate eastings and northings of a point layer
$x
$y

Calculate the length of a line layer
$length

Capitalise column values
upper(Field)
eg upper(Town)
Edinburgh becomes EDINBURGH

Camel case column values
title(Field)
EDINBURGH becomes Edinburgh
DUDDINGSTON LOCH becomes Duddingston Loch

Lower case column values
lower(Field)

Replacethis withthat in string
replace(string, replacethis, withthat)

Concatenate string a and string b
Concatenate a || b

Division and next line Multiplication
a/b
a*b

area/10,000 – divides area field by 10,000 (eg going from m2 to Hectares

Remove decimals from a field
toint(area)
eg 7954.235 becomes 7954 and 456525.325 becomes 456525

Index a set of polygons
$rownumber

Functions and Operators Official Notes for Field Calculator

Connecting to PostgreSQL 9.3 from QGIS 2.8.1 – local host

First ensure that you have both Postgres and QGIS installed on your machine.

In order for you to be able to connect to Postgres from QGIS on local host you must ensure 2 things. Firstly that the PostGIS plugin has been installed on your laptop AND secondly that you have included the postgis extension in each database that you wish to connect to. Without enabling the extension in the database you won’t be able to connect OR import shape files. Installation of PostGIS is often a default during the install of postgres but you can check whether this was completed correctly by using the Application Stack Builder, a small program that is installed with later versions of postgres.

I navigated to this on the win 8.1 machine I was using by using search.

Opening application stack builder you will be presented with the following.

ApplicationStackBuilder

Expand the spatial extensions tree to identify if you already have the PostGIS plugin installed – if not – select as appropriate the plugin and you will be prompted to install. You will need an internet connection for this. Above you can see that my plugin was already installed.

Next you need to add the PostGIS extension to each Postgres database you wish to link to from QGIS this is done through PG Admin.

This is something that both myself and a colleague got caught out by and it took me an hour of searching to find how to fix it.

Below I have a database called GISData which I have just created. You will note there is only one object within the expanded extensions tree. You will not be able to connect to a database that does not include PostGIS in its list of extensions
CreateExtension

Hi-light the database you want to spatially enable then go to Tools – Query Tool( Ctrl + E will do the same). In the above picture I’ve slightly jumped the gun. To add the extension to the database type.

CREATE EXTENSION postgis

Run the query by selecting the green right arrow
There will be a short delay and then upon refresh of the connection postgis should appear in the list of extensions.

CreateExtensionCreated

You can now close the Postgres administrator and return to QGIS where you should be able to setup the connection to the database.

Parameters should be similar to below and it is useful to test the connection prior to saving.

SettingupthePostGISconnection

Setting up a Blank SQL Server Spatially enabled Table using Microsoft SQL Server Management Studio 2008R2 Express and displaying it in QGIS 2.8.1

Programs used;

1- SQL Server 2008R2 Express
2- SQL Server Management Studio 2008R2 Express
3- QGIS

The example uses UK national grids coordinates to create a Triangle Polygon in a SQL Server Table

I’ve previously written that while we’ve had spatially enabled SQL Server for over 5 years I constantly come across line of business applications that although using SQL Server have not and do not intend to spatially enable the application. This is undoubtedly because of the difficulty in re-designing legacy systems actively in use and because the benefits although significant are not generally requested by all but the most knowledgable of colleagues.

While I understand this legacy system reasoning spatially enabled databases are the future so its just a matter of when and not if an application will require alteration. Understanding it in this context makes it really a requirement to start seriously planning for its inclusion.

Developerers creating new applications however should always consider spatially enabling relevant tables from the start even if it is not specked by the client/colleague. It being so much easier to spend a couple of minutes future proofing the schema of a new born database rather than hours trying to retrofit a live in production back end.

Firstly it’s important to understand what a geodatabase in SQL Server actually is.
Really it is a normal database which has one table that has a field that has a geometry or geography value type. In this example I will use desktop QGIS 2.8.1 to display the resulting geometry but any other digital mapping package that can link to SQL Server could be used. SQL Server also has a very rudimentary Mapping Display but you will need something better if you want to manipulate boundaries visually.

Many digital mapping products have plugins that will create Geodatabases and tables however I haven’t seen one for QGIS. I really wanted to be able to create spatial SQL tables on my own without recourse to paid tools directly in SQL Server Management Studio. So here’s my method of creating blank polygon table whose geometry is ready to be read and edited in QGIS or any other digital mapping system just using SQL Server Management Studio Express 08R2.

Steps
1. Create a new Table
2. Ensure the table has an identity Key that increments
3. Create a geometry column
4. Write a query that updates the geometry column

UPDATE T001Sites SET Coordinates=geometry::STGeomFromText(‘POLYGON((301804 675762,295789 663732,309581 664870,301804 675762))’,27700)

You will note that there are four coordinates here (each coordinate being a pair of numbers )
The first coordinate and last are the same this is required by SQL to ensure that the polygon is closed.

The 27700 number is the Spatial Reference System Identifier (SRID) – it is a unique value used to unambiguosly identify projecttion. 27700 is the unique identifier for the United Kingdom National Grid and represents the coordinates my example refer to. The spatial reference identification system is defined by the European Petroleum Survey Group (EPSG) standard which is a set of standards developmed for cartography surveying etc and owned by the Oil and Gas Producers Group list here; http://www.epsg-registry.org/

The above coordinates display a triangle in West Lothian near Edinburgh

5. Set up the connection to SQL Server Instance

Ensure the box marked “Only look in the geometry_columns metadata table” checkbox is unchecked. By default this is checked and if the geometry_columns table does not exist you will get an error message.

QGIS-SSMS-Connection

6. Display the table and edit as appropriate.

Select the table then hit the Add button

QGIS-SSMS-TableReadyforDisplay

And here is the SQL Server table in QGIS ready to be added to or edited.
QGISshowingSQLServerPolygon

Connecting to SQL Server – authentication and QGIS

Within QGIS when you set up a connection to a MS SQL Server instance you are presented with two options. Here’s a bit of clarification on what the two options entail.

* Trusted connection – this is the same thing as using Windows Authentication and authentication is managed by the domain and authorization is handled by SQL Server – This could be handled by an Active Directory Security Group.

* Login – SQL Server can also use its own logins such as a user. These are both authenticated and authorized by SQL Server. They are only viable if SQL Server is configured to run in Mixed Authentication mode.

image001

Introduction to Basic Printing QGIS 2.2

Creating maps that you can pass on to others is often a central and regular requirement if not in paper format then in a digital format that can be e-mailed or printed out. Here’s a quick reference for myself as much as anything else.

To get into the print composer you can create a completely new print composition or alternatively load an existing print composition – Generally the 5th icon in the main menu bar will take you there (should be a white landscape rectangle with a star will give you access to the print composer , demonstrated below;

The 6th icon can be used to get to an existing declared print composition.

PrintComposerDemonstrationEnvironment

Now in the first instance you are going to want to navigate around the map and ensure that the map you wish to produce has the correct extents. In the below image on a 2 screen image I show QGIS v2.2 open on the left screen and the print composer open in the right. To move the composition area around go to the map window within the main program and navigate accordingly. Then within the print composer window hit the command button titled

Set to map canvas extent

This will re-draw your composition with an interpreted boundary defined directly from you map window. You can enforce scale in the item properties. Similarly after changing layers you will need to ensure that you hit the above button again when you want the composition to reflect the layers within the map window.

QGIS – Free GREAT Digital Mapping Software

windglobe A map showing winds over the Atlantic

Looking for a desktop digital mapping package? You really need to check out QGIS it is an absolutely excellent open source geographical information system. At the time of writing the latest version was QGIS 2.4 – the below tips were taken from research into windows version of QGIS 2.2

Full program available here.
Link to www.qgis.org site (English)

Tip : Navigation – Magnification – Plus or Minus mangifier Icons or wheel scroll
Tip : Navigation – Scroll – cursor keys or alternatively the hand icon or hold down the space bar and movement of the mouse when pointer is in the map window.
Tip : Projection – CRS stands for Coordinate Referencing System – lots of different ways of showing what is essentially the surface of a sphere on a flat surface – and more generally referred to as map projection – you will remember from geography. For most UK maps the coordinates are often in Ordnance Survey UK Grid therefore you want the properties of Coordinate Referencing System of the project to be OSGB and you want the coordinate referencing system of the individual layers to be OSGB as well. Once this is done the scaling will be correct and so will the measurement tools.
Tip : View / Panel – allows you to switch on and off menus – very good and very powerful
Tip : Graphical Record selection – Icon in the middle of the toolbar that has a number of differing options – it’s a drop down that allows different things for selection.
Tip : Attribute Record selection – Icon in the middle of the toolbar that allows for table attribute selection. Shows the table and this can be sorted properly.
Tip : Deselect Records – can individually de-select using the keyboard alternatively you can also use the de-select icon in the middle of the top of the screen.
Tip : Browser – brilliant for navigating through the directory and seems a lot quicker than going through the pop up individual menus on the left – for me anyway – additionally you can add an additional browser layer and transfer things between directories. It is an excellent alternative to the file dialogue manager.
Tip : View / Decorations – You can add things like scale bar and copyright to the map window here – very intuitive and nice finishing touch to your projects.
Tip : Labelling – Make scale dependent – highlight the layer you are interested in and right click. Now select the Labels option and within the Size section change the drop down from points to map units.
Tip : Labelling – Threshold the labelling – right click on layer and then go to the Rendering section and select scale based visibilty and adjust accordingly.

Above interpreted from the QGIS manual see:
Link to PDF version of QGIS v2.2 manual