Copying Rasters in PostGIS

I ran into a process where I wanted to create copies of rasters in PostgreSQL. While seemingly a simple process this took me a bit of work to figure out.

For my workflow I had three rasters, which all have the same size, and I want to load them into the same PostGIS table with three raster geometry columns. I don’t think this will work for different sized rasters since the rid’s will not match.

Three rasters:
raster1
raster2
raster3

Which I want to copy into:
merged_raster

First to create the merged raster table:

CREATE TABLE merged_raster

(

  rid serial NOT NULL,

  raster1 raster,

  raster2 raster,

  raster3 raster

);

Then to add the rid’s. These are the id’s of the tiles that the raster was split into when loading. If your tile size is large enough then you may only have one.

INSERT INTO merged_raster(rid)

(SELECT rid FROM raster1);

Then copying the actual data is straighforward (this assumes the raster column in the raster1 datasets is called rast):

UPDATE merged_raster m

SET raster1 = r.rast

FROM raster1 r

WHERE r.rid = m.rid;



UPDATE merged_raster m

SET raster2 = r.rast

FROM raster2 r

WHERE r.rid = m.rid;



UPDATE merged_raster m

SET raster3 = r.rast

FROM raster3 r

WHERE r.rid = m.rid;

Now I still have an issue that QGIS will not load these layers. It will always load the initial raster column no matter what is chosen.

QGIS OpenStreetMap Scales

Save as a text file ending in .xml like qgis_scales.xml

These are the scales OpenStreetMap tiles are rendered in for 96 dpi, so the map will look sharp on most monitors. These are the scales for the zoom levels.

The xml file can then be loaded into the project from:

Project> Project Properties…> General> Project scales

<qgsScales version="1.0">
    <scale value="1:554678932"/>
    <scale value="1:277339466"/>
    <scale value="1:138669733"/>
    <scale value="1:69334866"/>
    <scale value="1:34667433"/>
    <scale value="1:17333716"/>
    <scale value="1:8666858"/>
    <scale value="1:4333429"/>
    <scale value="1:2166714"/>
    <scale value="1:1083357"/>
    <scale value="1:541678"/>
    <scale value="1:270839"/>
    <scale value="1:135419"/>
    <scale value="1:67709"/>
    <scale value="1:33854"/>
    <scale value="1:16927"/>
    <scale value="1:8463"/>
    <scale value="1:4231"/>
    <scale value="1:2115"/>
</qgsScales>

Example:

1,000,000 (QGIS default):

1,083,357 (OSM wiki):

1,155,584 (From: 3liz):

Scales from:
OSM wiki

Zoom scales.

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:

QGIS Select Within Plugin 0.4

Runs through each geometry that you want to select from and tests if the centroid or the point of surface central point falls within the selecting geometry. If the central point falls within the selecting polygon, it is usually mostly inside the selecting geometry.

Or selecting based on percentage overlap. By default dissolving the selecting features first.

Useful if you are splitting up a polygon layer based on a polygon coverage layer, as each polygon will only end up in one of the selecting polygon areas. Unlike with an Intersects query which would return the geometry twice if it is on the border one on each side, or a Within query, which would not return the geometry at all. Very useful if your selection and selecting geometries have similar boundaries.

Version 0.4 brings a couple of new features to the Select Within plugin.

Most importantly it now supports mostly within and percentage within selections.

Also introduced is Pole of Inaccessibility within selections. As implemented in the QGIS core with the Polylabel algorithm: https://github.com/mapbox/polylabel

Point within selection examples. Centroid (red), Point on Surface (green), and Pole of Inaccessibility (blue):

Examples:

Centroid within:

Point on Surface within:

Pole of Inaccessibility within (1.0 tolerance):

50% within:

75% within:

Code:
GitHub Code Base

Issues and bug tracker:
GitHub Issues

Original Post:
Centroid Within Selection in QGIS