Categories
All Animations Maps OSGEO PostGIS QGIS Tutorials

My Strava Activities from the Same Place at the Same Time Animated

All of my activities from Strava if they started from the same place at the same time.

Strava Activites Animated

Longest activity time wise is an 11 hour walk.

Distance wise is a 200km helicopter ride.

Includes a couple of 100km cycles. But mainly walks on the beach, 5 km runs, and 10km cycles.

You can see the patterns from a number of locations that I have lived in.

How to

  1. Download your activities from Strava.
  2. Parse them into a PostgreSQL/PostGIS database with Python: https://github.com/HeikkiVesanto/StravaGarminParser/blob/master/parser_points_to_db.py
  3. Process the points to time series in the database: https://gist.github.com/HeikkiVesanto/3fbd55cda45394d069773a34ea244e4b
  4. Create animation in QGIS using the atlas generator.
Categories
All Maps OSGEO QGIS Tutorials

Mapping Strava Data

In 2018 we started a running club at work.

I created a quick script to parse the data on Strava to a ShapeFile, which can be easily animated with QGIS.

The script only works with Garmin files, GPX, TCX, and FIT.

Source: https://github.com/HeikkiVesanto/Strava-Garmin-Parser

Example:

Categories
All Ireland OpenStreetMap OSGEO PostGIS QGIS Tutorials

Creating OpenStreetMap History Visualisations

I created a couple of OSM visualisations for my talk at the OSGeo Ireland conference.

See: History of OpenStreetMap in Ireland

These are pretty easy to make, but take a fair bit of time. I did mine for Ireland, but should work with any part of the world.

Required software:

  • PostgreSQL with PostGIS
  • Python
  • QGIS
  • osmium-tools

This is the trickiest part, installing osmium-tools: here.

Data:

An OSM full history export. The best source for these is GEOFABRIK.

For Ireland:

http://download.geofabrik.de/europe/ireland-and-northern-ireland.html

Due to GDPR, you will have to log in with an OSM id to download the full history extracts. User ID’s are personal data.

Process:

The workflow is pretty simple. Osmium-tools provides pretty easy API access to the history files, where you can provide a data, and it will extract what OSM was like at that date. We simply need to loop through the desired dates we want to extract, and pipe the results into a workflow that loads the data into PostgreSQL. The final step is simply rendering in QGIS using the time manager plugin.

Python Script:

Github GIST:

https://gist.github.com/HeikkiVesanto/f01ea54cca499a6a144d18cf8909c940

The tables in the database will be:

  • lines
  • multilinestrings
  • multipolygons
  • other_relations
  • points

Each feature will be tagged with the date it is associated with.

Visualisation:

To visualise the data in QGIS we use simply use the excellent time manager plugin, filtering on the load_date field and with a monthly interval.

Result:

Categories
All Animations Featured OSGEO QGIS Tutorials

Storm Harvey QGIS Geometry Generator

Storm Harvey produced some extremely high levels of rainfall. Some areas of Texas received over 50 inches of rain over 9 days.

The National Oceanic and Atmospheric Administration (NOAA) provided some really great real time datasets to map the progress of the storm.

Among these were:
Hourly Precipitation
and
Hurricane Path

From these we can produce a GIF of hourly precipitation:

Hourly precipitation.

And total precipitation:
Hurricane Harvey Total Precipitation

Particularly the hurricane path was possible to create in QGIS using the Atlas Generator, and the excellent new:ish geometry generator. This can be found as an option for any layers symbology, as one of the renderers.

For my map I had a non spatial table that drove my atlas. This was a log table of all of the hours of precipitation I had loaded into my database. So I looped through each entry and showed the corresponding points of hourly precipitation for the corresponding hour. I also had hurricane path data as points for every 6 hours. So I could use the geometry generator to interpolate points in between known points.

While the query ended up being pretty long it is pretty straightforward.

It only needs to be run when the hour being generated does not end with a 00, 06, 12, or 18, because those are the positions I already know.

For the rest I need to generate two points. One for the previously known point, and one for the next known point.

Then I would create a line between those two points, measure the line, and place a point on the line x times one sixth of the way for the start of the line depending on the hour from the last known point.

Overall I am very impressed and happy with the result. With a bit of data defined rotation the storm progress looks great.

line_interpolate_point(
Make_line(
geometry(
case when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12', '18') then
    get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') )
else
    get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - (attribute(  @atlas_feature , 'id') % 100 % 6  ))
end)
,
geometry(
case
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 6 )
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('18') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 100 - 18 )
    when to_int(right(to_string(attribute(@atlas_feature , 'id')),2)) > 18 then
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 100 - 18)
    else
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 6)
end)
),
length(Make_line(geometry(
case when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12', '18') then get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') )
else get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - (attribute(  @atlas_feature , 'id') % 100 % 6  ))
end
)
,
geometry(
case
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('00', '06', '12') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 6 )
    when right(to_string(attribute(@atlas_feature , 'id')),2) IN ('18') then
        get_feature(  @layer_name , 'dtg', attribute(  @atlas_feature , 'id') + 100 - 18 )
    when to_int(right(to_string(attribute(@atlas_feature , 'id')),2)) > 18 then
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 100 - 18)
    else
        get_feature(  @layer_name , 'dtg',  attribute(  @atlas_feature , 'id') - ((attribute(  @atlas_feature , 'id') % 100 % 6)  ) + 6)
end)))
*
((attribute(  @atlas_feature , 'id') % 100 % 6) * 0.16666666666666666))

Categories
All OSGEO QGIS Tutorials

Mapping Google Location Data

A cool python script has been created that allows you to easily convert your google location (Takeout) data into a shapefile.

You can get your data from: Google Takeout
And you only need the “Location History – JSON format”

The conversion python script can be downloaded from: GitHub

The python script requires GDAL and its python bindings, but can be easily run if you installed QGIS using the OSGeo4W installer. From the advanced installer, under the Lib section.

instruct

Then using the OSGeo4W Shell.
shell

Run the command:

python "C:\FullPath_to_Python_Script\read_location_data.py" "C:\FullPath_to_Input_File\LocationHistory.json" "C:\output_path" output_file_name ESRI_Shapefile

Example:

python "C:\FilePath\Takeout\Location History\read_location_data.py" "C:\FilePath\Takeout\Location History\LocationHistory.json" "C:\FilePath\Takeout\Location History" output ESRI_Shapefile

Then just style it in QGIS as desired.
GoogleTakeOut

Categories
All OSGEO Tutorials

GIS to CAD using ogr2ogr – Part 3 – Point Annotation to Text in CAD

GIS to CAD using ogr2ogr – Part 1 – Shp to DXF with Contour Data
GIS to CAD using ogr2ogr – Part 2 – GML to DXF with OS MasterMap

MasterMap Topo Sample Data:

https://www.ordnancesurvey.co.uk/business-and-government/licensing/sample-data/discover-data.html

OS MasterMap has an annotation layer, which is simple to symbolise in a GIS program. But becomes more difficult in CAD software.

With ogr2org, when writing a DXF file, if you have an input point geometry, which has an OGR_STYLE attribute, it will be written as a text geometry when opened in CAD.

So for our MasterMap data we have one layer we want to convert to text:
CartographicText

ogrinfo os-mastermap-topography-layer-sample-data.gml CartographicText
OGRFeature(CartographicText):855
fid (String) = osgb1000000729425681
featureCode (Integer) = 10026
version (Integer) = 1
versionDate (String) = 2001-11-11
theme (StringList) = (1:Buildings)
changeDate (StringList) = (1:1999-09-08)
reasonForChange (StringList) = (1:Modified)
descriptiveGroup (StringList) = (1:Buildings Or Structure)
physicalLevel (Integer) = 50
anchorPosition (Integer) = 0
font (Integer) = 2
height (Real) = 1.5
orientation (Integer) = 2890
textString (String) = 5
descriptiveTerm (String) = (null)
make (String) = Manmade
POINT (290875.38 92635.81)

So for this we are primarily interested in “textString” and potentially “orientation”.

Lets see the layer as points first as a baseline:

ogr2ogr -f DXF CartographicText.dxf os-mastermap-topography-layer-sample-data.gml CartographicText

Screenshot[10]

Zoomed in:

Screenshot[11]

But lets try that as text. We will keep this simple and only take into account orientation and to a small extent height. Lets look at orientation:

Orientation – The orientation of text or symbol features for cartographic placement. This is measured in tenths of a degree anticlockwise from due east (0–3599).

So conversion to degree will be simple. Orientation/10

We can also take into consideration height as a multiplier.

And “textString” stores the text itself.

The command:

ogr2ogr.exe -f DXF CartographicText2.dxf os-mastermap-topography-layer-sample-data.gml CartographicText -sql "SELECT 'LABEL(f:""Arial"",s:""'||(height*800)||'"",t:""'||textString||'"",a:""'||(orientation/10)||'"",p:5)' AS OGR_STYLE, * FROM CartographicText" -dialect SQLITE

Full extent:

Screenshot[12]

Zoomed in:

Screenshot[13]

Explained:

Since this is run in windows, through the regular console, the escape character for quotes is two quotes “”‘. So a combination on ‘ ” and “”‘ we can accommodate all the required quotes.

f:””Arial””,s:””‘||(height*800)||'””,t:””‘||textString||'””,a:””‘||(orientation/10)||'””,p:5

f:””Arial””
Font: Arial

s:””‘||(height*800)||'””
Size: Multiplier of the height field, I am not sure what the units are, comments appreciated.

t:””‘||textString||'””
Text: textString column

a:””‘||(orientation/10)||'””
Align: In degrees

p:5
Position: the OS position and the ogr2ogr style position are slightly different, so better placement could be achieved with some pre-processing

GIS to CAD using ogr2ogr – Part 1 – Shp to DXF with Contour Data
GIS to CAD using ogr2ogr – Part 2 – GML to DXF with OS MasterMap

Categories
All OSGEO QGIS Tutorials

GIS to CAD using ogr2ogr – Part 2 – GML to DXF with OS MasterMap

GIS to CAD using ogr2ogr – Part 1 – Shp to DXF with Contour Data
GIS to CAD using ogr2ogr – Part 3 – Point Annotation to Text in CAD

For this example we are using Ordnance Survey MasterMap Topology Layer data.

MasterMap Topo Sample Data:

https://www.ordnancesurvey.co.uk/business-and-government/licensing/sample-data/discover-data.html

Now we know that we can maintain an attribute through layers, as we saw in the shp to DXF example, the export of MasterMap should be straightforward.

Let’s first see what the GML file contains.

ogrinfo -so os-mastermap-topography-layer-sample-data.gml
Had to open data source read-only.
INFO: Open of 'os-mastermap-topography-layer-sample-data.gml'
using driver 'GML' successful.
1: TopographicArea (Polygon)
2: CartographicText (Point)
3: CartographicSymbol (Point)
4: BoundaryLine (Line String)
5: TopographicPoint (Point)
6: TopographicLine (Line String)

So we have 6 layers in total.

For MasterMap in CAD we will be mainly interested in CartographicText, TopographicPoint, and TopographicLine.

Lets start with TopographicLine.

ogrinfo -so os-mastermap-topography-layer-sample-data.gml TopographicLine

Nothing too useful.

A bit more details:

ogrinfo os-mastermap-topography-layer-sample-data.gml TopographicLine

Sample:

OGRFeature(TopographicLine):185
fid (String) = osgb1000000347615024
featureCode (Integer) = 10019
version (Integer) = 5
versionDate (String) = 2005-03-30
theme (StringList) = (2:Buildings,Land)
accuracyOfPosition (String) = 1.0m
changeDate (StringList) = (4:1994-01-26,2003-11-10,2004-02-19,2005-01-05)
reasonForChange (StringList) = (4:Modified,Attributes,Attributes,Attributes)
descriptiveGroup (String) = Building
physicalLevel (Integer) = 50
physicalPresence (String) = Obstructing
descriptiveTerm (String) = Outline
make (String) = Manmade
nonBoundingLine (String) = (null)
LINESTRING (...)

For this feature the “descriptiveGroup”” seems the most useful, and from reading the os-mastermap-topography-layer-user-guide.pdf the best would be either a combination of descriptiveGroup and descriptiveTerm or using the featureCode. Since this is a simple conversion we will just use a combo of descriptiveGroup and descriptiveTerm to create our DXF layers.

I will be using || for concatenation, which works with the SQlite SQL dialect.

ogr2ogr -f DXF TopographicLine.dxf os-mastermap-topography-layer-sample-data.gml TopographicLine -sql "select descriptiveGroup || ' - ' || descriptiveTerm as Layer from TopographicLine" -dialect SQlite

Result:

layer names ignored in combination with -sql.
ERROR 1: No known way to write feature with geometry 'None'.
ERROR 1: Unable to write feature 0 from layer SELECT.

ERROR 1: Terminating translation prematurely after failed
translation from sql statement.

Not quite. Seems to be missing geometry, perhaps a SQL select issue.

This can be tested with:

ogrinfo os-mastermap-topography-layer-sample-data.gml TopographicLine -sql "select descriptiveGroup || ' - ' || descriptiveTerm as Layer from TopographicLine" -dialect SQLITE

Result:

OGRFeature(SELECT):14634
Layer (String) = Building - Outline

OGRFeature(SELECT):14635
Layer (String) = Building - Outline

So we do not have any geometry. Lets bring that in.

ogr2ogr -f DXF TopographicLine.dxf os-mastermap-topography-layer-sample-data.gml TopographicLine -sql "select descriptiveGroup || ' - ' || descriptiveTerm as Layer, * from TopographicLine" -dialect SQLITE

Geometry looks good:

Screenshot[6]

But if we check the attributes in QGIS:

Screenshot[7]

We can see that all of the attributes that are not 0 have both a descriptiveGroup and a descriptiveTerm, which was not what we can see in the ogrinfo summary. So our SQL statement is cutting some out.

Try again:

ogr2ogr -f DXF TopographicLine2.dxf os-mastermap-topography-layer-sample-data.gml TopographicLine -sql "select descriptiveGroup ||' - '|| coalesce(descriptiveTerm,'') as Layer, * from TopographicLine" -dialect SQLITE

Looking better:

Screenshot[8]

But it won’t open in AutoCAD DWG TrueView. Lets try running it through a ShapeFile format first before the DXF conversion.

ogr2ogr TopographicLine.shp os-mastermap-topography-layer-sample-data.gml TopographicLine -sql "select descriptiveGroup || ' - ' || coalesce(descriptiveTerm,'') as Layer, * from TopographicLine" -dialect SQLITE

ogr2ogr -f DXF TopographicLine3.dxf TopographicLine.shp

Success:

Screenshot[9]

No indication of why a direct GML to DXF conversion would hang TrueView, and your mileage with other CAD software may vary. But ShapeFile is a very simplified geometry format, so perhaps running through that helps with some more complex geometry in the GML. Hard to say with no errors from TrueView, just a stuck program.

Repeat for point:

ogr2ogr -f DXF TopographicPoint.dxf TopographicPoint.shp

ogr2ogr TopographicPoint.shp os-mastermap-topography-layer-sample-data.gml TopographicPoint -sql "select descriptiveGroup || ' - ' || coalesce(descriptiveTerm,'') as Layer, * from TopographicPoint" -dialect SQLITE

GIS to CAD using ogr2ogr – Part 1 – Shp to DXF with Contour Data
GIS to CAD using ogr2ogr – Part 3 – Point Annotation to Text in CAD

Categories
All OSGEO QGIS Tutorials

GIS to CAD using ogr2ogr – Part 1 – Shp to DXF with Contour Data

GIS to CAD using ogr2ogr – Part 2 – GML to DXF with OS MasterMap
GIS to CAD using ogr2ogr – Part 3 – Point Annotation to Text in CAD

The power of GDAL, and specifically ogr2ogr is pretty impressive. This conversion is from shp to DXF, which is a somewhat universal CAD format so further conversion should be possible.

This post will cover contour export while maintaining 3D elevation, in addition to contour values as layers in CAD. The data used is OS terrain 50.

OS Terrain 50:

https://www.ordnancesurvey.co.uk/opendatadownload/products.html

In QGIS:

Screenshot[1]

Contours in 3D:

ogr2ogr -f DXF contour_zfield.dxf SX99SW_line.shp -zfield PROP_VALUE

With the -zfield creating the 3d elevation.

Great result:

Screenshot[3]

The alternative is to just store the z-value as layers.

ogr2ogr -f DXF contour_layer.dxf SX99SW_line.shp -sql "SELECT PROP_VALUE AS Layer FROM SX99SW_line"

Layers work great:

Screenshot[2]

With the ogr2ogr DXF driver, if you have an input column called “Layer” then it will be used to group features as a layer in DXF. We use a SQL query to achive this. Prop_Value is the height field in my input data.

And putting them all together:

ogr2ogr -f DXF contour_zfield_layer.dxf SX99SW_line.shp -zfield PROP_VALUE -sql "SELECT PROP_VALUE AS Layer FROM SX99SW_line"

Result not as expected, flat output:

Screenshot[4]

Adding our SQL select statement removes our zfield attribute as such ogr2ogr cannot access it. Lets resolve this:

ogr2ogr -f DXF contour_zfield_2_layer.dxf SX99SW_line.shp -zfield PROP_VALUE -sql "SELECT PROP_VALUE AS Layer, * FROM SX99SW_line"

Excellent:

Screenshot[5]

Layers and height.

GIS to CAD using ogr2ogr – Part 2 – GML to DXF with OS MasterMap
GIS to CAD using ogr2ogr – Part 3 – Point Annotation to Text in CAD

Categories
All Tutorials

Setting up your GDAL and OGR Environmental Variables

Table of Contents

    Windows 11

    Open the start menu, and search for: environment variables

    Choose the Edit the system environment variables

    Choose Environment Variables.

    Edit the Path Variable by adding an entry for you path to ogr2ogr.

    I would recommend installing ogr2ogr using the OSGeo4W network installer and the default path. If you have done this the path is: C:\OSGeo4W\bin

    You also need to add 2 new entries into the base environmental variables, one for:

    GDAL_DATA

    With path:

    C:\OSGeo4W\share\gdal

    One for:

    PROJ_LIB

    With path:

    C:\OSGeo4W64\share\proj

    Adjust the paths accordingly if you have not used the default install locations.

    Windows 10

    I recently upgraded to Windows 10, and have been setting up my GIS workspace settings.

    Setting up your Windows path environmental variables can easily be done through the command line.

    These commands can be run directly in the command line and don’t require a computer restart to take effect. To run the command prompt, simply open up the windows search (Windows Key + S) and search for “cmd”:

    Screenshot[3]

    Then run the command to add the path.

    Update: PROJ_LIB also needs to be set, to prevent “PROJ: proj_create_from_wkt: Cannot find proj.db” errors.

    Determine where your QGIS is installed. I recommend using the OSGeo4W network installer and using the default location. Which would be: C:\OSGeo4W\bin

    If it is installed there, just run:

    setx PATH "%PATH%;C:\OSGeo4W\bin"
    setx GDAL_DATA "C:\OSGeo4W\share\gdal"
    setx PROJ_LIB "C:\OSGeo4W\share\proj"

    In the past a 64 bit install was installed by default into: C:\OSGeo4W64\bin

    If it is installed there run:

    setx PATH "%PATH%;C:\OSGeo4W64\bin"
    setx GDAL_DATA "C:\OSGeo4W64\share\gdal"
    setx PROJ_LIB "C:\OSGeo4W64\share\proj"
    Screenshot[6]

    Restart your command prompt.

    Screenshot[7]

    Success.

    No more:
    ‘org2ogr’ is not recognized as an internal or external command, operable program or batch file.

    Categories
    All OSGEO QGIS Tutorials

    Georeferencing vector data using QGIS and ogr2ogr

    Update 2022

    This can also be achieved in QGIS now.

    Since QGIS 3.26 this can now be achieved with the georeferencer, same tool that is used of raster georeferencing:

    QGIS now supports georeferencing vector layers in the georeferencer tool. This allows vector layers without spatial referencing to be interactively georeferenced, or layers with referencing to be re-referenced, in a similar manner to raster data. Georeferencing occurs in a task, so QGIS should remain responsive, even with large datasets.

    Changelog

    This post will cover:

    1. Digitising rasters using QGIS.
    2. Loading data to PostGIS using SPIT.
    3. Reprojecting raster data using QGIS.
    4. Lines within a polygon selection in PostGIS.
    5. Georeferencing vector files using ogr2ogr, search for “The georeferencing:” to skip the other steps.

    Georeferencing vector data has long been very difficult using the open-source GIS stack. That has all changed with the release of GDAL 1.10.0. Now georeferencing can be done using the vector translator ogr2ogr.

    About half a decade ago I wrote my undergraduate dissertation on: “The influence of the Blaeu Atlas of Scotland on subsequent maps of Scotland”. While I was very proud of the work at the time, and the grade was good, the markers comments could be paraphrased as: A very good literature review on the Blaeu Atlas, however somewhat weak on analysis”. After a Master’s degree I could not agree more. At the time I lacked the knowledge and experience of working with projections. The release of GDAL 1.10.0 gives the perfect opportunity to return and correct those mistakes.

    One of my favourite aspects of the Blaeu analysis was a comparison of the Blaeu 1654 Atlas of Scotland coastline to the coastline of subsequent maps of Scotland. I originally did this in GIMP and manually resized the outlines to close approximations. The end result looked good, but it lack real scientific rigour.

    Original comparison:

    1654_Blaeu_map_of_scotland_outline

    A huge thanks to the National Library of Scotland for their work in digitising the Blaeu Atlas of Scotland. A great online version and overview can be found at:  http://maps.nls.uk/atlas/blaeu/

    With the outline used today found at: http://maps.nls.uk/view/00000383

    Originally the maps were provided as .jpeg files with a simple viewer used for zooming. This meant that editing the URL would allow one to retrieve a very large version of the of the map directly. Now the maps are served using a a very nice javascript map viewer and digital copies can be purchased.

    First we need to create a vector outline from the raster image. This is a simple case of adding the .jpeg image in as a raster file. When prompted for a CRS I chose EPSG:4030 “Unknown datum based upon the WGS 84 ellipsoid” this is just used for meta-data at this point. What you choose does not matter. We just want a vector file in cartesian co-ordinates. We create a new shapefile layer and trace outline. My co-ordinate capture will be very rough. If you are working with more modern maps, you should be as careful as required:

    Screenshot from 2014-01-20 20:43:20

     

    Whenever I needed a break I would end the line. Then I could start a new line snapped to the end of the last one. In the end you can select all of the line segments and “merge selected features” from the advanced digitising toolbar. This would likely have been enough, but I was concerned there would be two overlapping nodes and that I would have a multi-part line. So I exploded the lines into points using the Vector>Geometry Tools>Extract nodes tool. Then Using the handy Points2One plugin (thanks Pavol Kapusta) I could stitch the nodes back together into a single line. If you’re going to do something might as well do it correctly. In the end we have 5116 nodes:

    Screenshot from 2014-01-20 23:45:09

    With the Blaeu outline captured we need a basemap to georeference it to. In the end we want to re-projectiong it into the WGS 84 / World Mercator projection (EPSG:3395). Mercator will provide a good base of what the original mapper would have been trying to capture. In the UK we can turn to the Ordnance Survey for a basemap, however that basemap needs to be reprojected to EPSG:3395. This is simply done in within QGIS by selecting the “Enable on-the-fly CRS transformation” to EPSG:3395 from Project>project Properties>CRS tab.

    From the Ordnance Survey (OS) we have a few open datasets to choose from:

    • OS BoundaryLine (High Water Polyline)
    • OS BoundaryLine (european electoral regions)

    We can see how they line up with the raster products:

    • OS MiniScale
    • OS 250K

    We can see that the High Water Polyline follows MiniScale much better than the electoral regions:

    Screenshot from 2014-01-21 20:13:23
    Screenshot from 2014-01-21 20:13:41

    It doesn’t line up perfectly with the 250k raster (250k is probably too accurate for this study):

    Screenshot from 2014-01-21 20:19:07

    But on the whole it will serve for this purpose, in addition Scotland can be easily extracted:

    Screenshot from 2014-01-21 21:12:47
    Screenshot from 2014-01-21 21:04:42

    Now we have the High Water Polyline, which looks good, but we only want mainland Scotland. Easy right?

    Well… slighlty more complex. Unfortunately even the Scottish high water lines consisted of 16030 line segments. Select by location in QGIS using the European Electoral Regions (polygon) took too long (hours rather than minutes). So I had to resort to PostGIS to do the selection.

    SPIT to PostGIS:

    Screenshot from 2014-01-21 21:51:16

    For setting up a PostGIS database see my previous post: Setting up PostgreSQL and PostGIS on Linux Mint.

    The query in PostGIS/PostgreSQL to select all of the lines within a polygon (we use the OS boundary line electoral regions):

    CREATE TABLE scotland_selection AS
    SELECT os_high_water_mark.* FROM os_high_water_mark,european_scotland
    WHERE ST_Intersects(european_scotland.the_geom,os_high_water_mark.the_geom);

    Result:

    Query returned successfully: 4994 rows affected, 764779 ms execution time.

    The 12 minute run time says a lot about my poor computer.

    Then the lines are joined together, the geometries simplified, and multipart converted to singlepart. This allows us to get rid of the final islands. With the result back in QGIS:

    Screenshot from 2014-01-22 18:03:40

    Closeup:

    Screenshot from 2014-01-23 00:04:17

    Final result projected into WGS 84 / World Mercator projection (EPSG:3395), which really shows the distortion caused by the British National Grid (EPSG:27700):

    Screenshot from 2014-01-23 00:12:45

    We now have a cartesian vector file ready for georeferencing, an OS raster background map to use for georeferencing (OS MiniScale) and an OS vector outline BoundaryLine (High Water Polyline) to use for the final comparison.

    The georeferencing:

    We open up two QGIS projects. In one we have our cartesian vector file. In the other we have our basemap:

    Screenshot from 2014-01-23 00:16:01

    Using the co-ordinate capture tool, we capture the same point from both projects. Basically we are just capturing co-ordinate pairs from the two maps. The “Copy to clipboard” feature comes in very handy here.

    Screenshot from 2014-01-21 18:32:34

    These captured co-ordinates should be pasted into a text file. Something like this:

    Dumfries and Galloway
    1571.29165 -4217.06023 1571.29165 -4217.06023
    -551760.015 7302342.218 -551760.015 7302342.218
    
    Scotland England border
    2729.71384 -4421.27506 2729.71384 -4421.27506
    -338736.496 7324447.924 -338736.496 7324447.924
    
    Kintyre:
    1477.31668 -3328.36427 1477.31668 -3328.36427
    -603897.130 7496225.431 -603897.130 7496225.431

    Of note here is that the co-ordinate capture tool captured the co-ordinate twice for me. How it works is that it captures one in the original CRS and the selected one. It is also important to keep your co-ordinates in order. I captured the cartesian co-ordinates on the first row and the WGS84 ones on the second row. Also I have attached memorable names to the points, in case there is an issue, the culprit point can be easily be identified.

    For my project 23 points were captured:

    Screenshot from 2014-01-23 00:48:57

    After the points are captured we switch to ogr2ogr to perform the actual reprojection and to convert the OS outline into World Mercator.

    Scotland Mainland Outline from Boundary Line High Water Mark to World Mercator (while I’m converting I will convert it to a GeoJSON for viewing in leaflet (this will be covered by a later post) and strip out any attributes to decrease the file size):

    ogr2ogr -progress -f GeoJSON blaeu_outline.geojson -a_srs EPSG:3395 -dsco ATTRIBUTES_SKIP=YES -order 1 -gcp 1571.29165 -4217.06023 -551760.015 7302342.218 -gcp 2729.71384 -4421.27506 -338736.496 7324447.924 -gcp 1477.31668 -3328.36427 -603897.130 7496225.431 -gcp 2011.25239 -3353.35609 -536774.862 7514805.007 -gcp 2021.23961 -3783.80519 -516238.228 7418369.368 -gcp 1386.05253 -2433.53330 -632410.002 7624335.089 -gcp 1609.76622 -1681.49577 -635963.491 7783695.427 -gcp 1732.60900 -1079.26651 -599335.214 7983647.549 -gcp 2102.13607 -655.80846 -556830.013 8062917.700 -gcp 2906.10714 -611.86470 -339110.443 8067291.226 -gcp 2974.02022 -753.68320 -340067.151 8031482.985 -gcp 2858.16849 -1518.70411 -423164.137 7901097.254 -gcp 2617.47653 -1878.74333 -470179.537 7823193.830 -gcp 3788.47786 -1882.73821 -224168.724 7866382.395 -gcp 3920.30914 -2064.50559 -200661.024 7828114.046 -gcp 3271.13996 -3004.30282 -288951.571 7635405.575 -gcp 3216.21026 -3254.98199 -289771.607 7578549.743 -gcp 2608.98740 -3372.33181 -415578.804 7531397.670 -gcp 3120.33297 -3431.75575 -312337.131 7534907.369 -gcp 3433.93162 -3683.43365 -238109.336 7506113.225 -gcp 3495.85237 -3831.74384 -222391.979 7477821.982 -gcp 2050.70190 -4455.44562 -457276.805 7286175.723 1654_Blaeu_map_of_scotland.shp

    The key here is the -gcp tag. For each control point (co-ordinate pair) we have one tag:

    Dumfries and Galloway
    1571.29165 -4217.06023 1571.29165 -4217.06023
    -551760.015 7302342.218 -551760.015 7302342.218

    Turns into:

    -gcp 1571.29165 -4217.06023 -551760.015 7302342.218

    So: -gcp ungeoref_x ungeoref_y georef_x georef_y

    One -gcp tag for each point.

    -order 1

    Means that it will use the a first order polynomial transformation, so affine. For myself the less distortion the better.

    The end result is great, it really shows how advanced the Blaeu Atlas was for a 1654 map:

    blaeu_os_comparison

    If you are interested in buying a historically amazing map:

    Blaeu Atlas Maior of 1665 – Including the atlas of Scotland.

    Also please visit the NLS, without who this would have been impossible, they also have great changing exhibits:

    National Library of Scotland

    NLS Exhibitions