Robin's Blog

Some easy recipes

This is a bit different to what I usually post on this blog (which is usually technical content about GIS, remote sensing, Python and data analysis – see a few example posts), and this is part of April Cools – a group of bloggers writing posts that are unusual for them, but still good content.

So, on that note, here are a few fairly simple recipes that we often make in our household. None of them are particularly fancy, but they all taste nice – and my family like them.

Sausage, apple and pepper thing – in the process of being cooked

Cheat’s Carbonara

Serves approx 3 people.

  1. Chop chestnut mushrooms (half a tub) and bacon (approx. 4-6 rashers).
  2. Fry bacon in a little oil, add mushrooms after a couple of minutes
  3. Meanwhile, cook tagliatelle pasta in a separate saucepan, according to packet instructions
  4. After about 5 minutes, add double-cream and a good sprinkling of pepper to the frying pan, and cook until the cream thickens and absorbs some of the bacon fat.
  5. Mix in the cooked pasta, and serve with cherry tomatoes.

Cheat’s Curry

Serves approx 4 people.

  1. Dice some chicken breast (or any other meat!), and fry until lightly browned
  2. Add a tablespoon of curry paste (Korma paste or Thai Green Curry paste work well, but I imagine others would work too). Stir to coat.
  3. Add a tin of coconut milk, stir to mix in the coconut solids. Add seasonings (mixed herbs, pepper, salt) to taste.
  4. Cook until sauce thickens, and serve with rice
  5. Optional: Add diced peppers or other vegetables to the pan at the start.

Coronation Rice

(lovely as a side for a picnic or barbeque)

  1. Put a heaped spoonful of turmeric into a saucepan full of water, and cook rice in the saucepan according to the packet instructions.
  2. Drain rice and cool under running cold water
  3. Cut up dried apricots into very small pieces, and mix these and a few handfuls of sultanas into the rice
  4. Add curry power, mayonnaise and mango chutney to the bowl and stir until well-mixed (you may find it easier to mix these three ingredients separately and then add to the main bowl once combined)
  5. Refrigerate and leave to develop a deeper flavour.

Sausage, Apple and Pepper Thing

(Imaginative name, isn’t it!)
Serves approx 4 people.

  1. Cut 3 Braeburn apples into small slices. Put in a frying pan with a little oil and fry until they’re slightly golden in colour.
  2. Add chopped peppers and cook briefly.
  3. Sprinkle sage into the pan to taste, and add salt and pepper as required.
  4. Chop sausages into 2cm long chunks and add to the frying pan. Cook until lightly browned.
  5. Add apple juice to the pan, until it comes about half-way up the contents. Crumble in a vegetable or chicken stock cube.
  6. Cook until sauce has reduced significantly, adding more liquid if needed, and ensuring that the sausages are fully cooked.
  7. Serve with rice.

Explore UK place names on a map – my new web app

As an Easter present for you all, I’ve got a new web app you can play with.

It lets you search for UK place names – things like ‘ends with burgh’ or ‘starts with great’ or ‘contains sea’ – and then plot them on an interactive map. Once you’ve done that, you can even share links to your settings so other people can see your findings.

Have a look at the web app now, and please tweet/toot/email me with any interesting things you find.

A few examples:

The data comes from the Ordnance Survey Open Names dataset, and it was processed on my local machine using Python to create some filtered data that could be used in the webapp, which all runs client-side.


My geospatial PDF talk at FOSS4G 2021

This is only about 3 years late – but I gave a talk at FOSS4G 2021 on geospatial PDFs. The full title was:

From static PDFs to interactive, geospatial PDFs, or, ‘I never knew PDFs could do that!’

The video is below:

In the talk I cover what a geospatial PDF is, how to export as a geospatial PDF from QGIS, how to import that PDF again to extract the geospatial data from it, how to create geospatial PDFs using GDAL (including styling vector data) – and then take things to the nth degree by showing a fully interactive geospatial PDF, providing a UI within the PDF file. Some people attending the talk described it as "the best talk of the conference"!

A few relevant resources are below:


Fixing SQL error ‘invalid reference to FROM-clause entry for table’

Here’s a bit of SQL I wrote recently that had an error in it which I struggled to find. The error I got was ERROR: invalid reference to FROM-clause entry for table "roads" and I did some Googling but nothing really seemed to make sense as a cause of this error.

My SQL looked like this:

WITH args AS (
    SELECT
        ST_TileEnvelope(z, x, y) AS bounds,
        ST_Transform(ST_MakeEnvelope(l, b, r, t, 4326), 27700) AS area
),
mvtgeom AS (
    SELECT
        ST_AsMVTGeom(
            ST_Transform(ST_Intersection(roads.geom, args.area), 3857),
            args.bounds
        ) AS geom
    FROM
        roads, args
    JOIN floods ON ST_Intersects(roads.geom, floods.wkb_geometry)
    WHERE
        roads.geom && args.area
        AND roads.geom && args.bounds
)
SELECT
    ST_AsMVT(mvtgeom, 'default') INTO result
FROM
    mvtgeom;

This looks a bit weird, but it’s a slightly re-written version of the SQL inside a function that was written to produce Mapbox Vector Tiles (MVT) for use in an online webmap. The function is picked up by pg_tileserv and called appropriately when tiles are needed.

What I’m doing in the code is setting up a couple of parameters – the tile envelope and the selected area, and then running a fairly simple join with extra WHERE conditions, converting the outputs to MVT geometries and then the whole thing to a MVT itself.

The problem is inside the mvtgeom AS section – let’s simplify the SQL there a bit so it is easier to look at:

SELECT
        ST_Intersection(roads.geom, args.area), 3857) AS geom
FROM
        roads, args
        JOIN floods ON ST_Intersects(roads.geom, floods.wkb_geometry)
WHERE
        roads.geom && args.area
        AND roads.geom && args.bounds

I’ve got rid of all the MVT related stuff, so we have a more standard PostGIS query: we’re looking for the intersection of the roads and the area, where the roads intersect the floods, and the roads are inside the area and inside the tile boundaries (both defined in the args CTE). Possibly the error is a bit clearer now – maybe…

The answer is, it’s all about the placement of the JOIN clause. I’d never really properly thought about this (even though it’s a bit obvious!), but the JOIN statement refers to the table listed immediately before it. Therefore, if we reformat the above slightly we’re selecting:

FROM
roads,
args JOIN floods ON ST_Intersects(roads.geom, floods.wkb_geometry)

That is, we are running the JOIN on the args table, not on the roads table. The slightly confusing error message (ERROR: invalid reference to FROM-clause entry for table "roads" ) means that this JOIN clause makes no sense when run against the args table, as there’s no such thing as roads in that context.

All we need to do to fix this is move the JOIN clause before args:

FROM
roads JOIN floods ON ST_Intersects(roads.geom, floods.wkb_geometry),
args

Remember the comma after the JOIN clause – we’re still in the list of tables/things that we’re selecting FROM.

In general, the format of a FROM clause is:

FROM
<tablename> [JOIN clause],
<tablename>,
…

For context, I was running this SQL in Postgres, with the PostGIS extension. I’m not an expert in SQL, so please let me know if I’ve made any errors in this post.


One reason for getting a ‘No HTTP triggers found’ error when using Azure Functions with Python V2 programming model

Summary: It might be because there is an exception raised when importing your function_app.py – for example, caused by one of your import statements raising an exception, or a parsing error caused by a syntax error.

I was deploying a FastAPI app to Azure Functions recently. Azure Functions is the equivalent of AWS Lambda – it provides a way to run serverless functions.

Since I’d last used Azure Functions, Microsoft have introduced the Azure Functions Python V2 programming model which makes it easier and cleaner to do a number of common tasks, such as hooking up a FastAPI app to run on Functions.

However, it also led to an error that I hadn’t seen before, and that I couldn’t find documented very well online.

Specifically, I was getting an error at the end of my function deployment saying No HTTP triggers found. I was confused by this because I had followed the documented pattern for setting up a FastAPI app. For reference, my function_app.py file looked a bit like this:

import azure.functions as func
from complex_fastapi_app import app

app = func.AsgiFunctionApp(app=app,
            http_auth_level=func.AuthLevel.ANONYMOUS) 

This was exactly as documented. But I kept getting this error – why?

I replaced the import of my complex_fastapi_app with a basic FastAPI app defined in function_app.py, this time copied directly from the documentation:

import azure.functions as func 
from fastapi import FastAPI, Request, Response 

fast_app = FastAPI()

@fast_app.get("/return_http_no_body") 
async def return_http_no_body(): 
    return Response(content="", media_type="text/plain") 

app = func.AsgiFunctionApp(app=fast_app, 
                           http_auth_level=func.AuthLevel.ANONYMOUS) 

Everything worked fine now and I didn’t get the error.

After a lot of debugging, it turns out that if there is an exception raised when importing your function_app.py file then Functions won’t be able to establish what HTTP triggers you have, and will give this error.

In this case, I was getting an exception raised when I imported my complex_fastapi_app, and that stopped the whole file being processed. Unfortunately I couldn’t find anywhere that this error was actually being reported to me – I must admit that I find Azure logging/error reporting systems very opaque. I assume it would have been reported somewhere – if anyone reading this can point me to the right place then that’d be great!

I’m sure there are many other reasons that this error can occur, but this was one that I hadn’t found documented online – so hopefully this can be useful to someone.


How to speed up appending to PostGIS tables with ogr2ogr

Summary: If appending to a PostGIS table with GDAL/OGR is taking a long time, try setting the PG_USE_COPY config option to YES (eg. adding --config PG_USE_COPY YES to your command line). This should speed it up, but beware that if there are concurrent writes to your table at the same time as OGR is accessing it then there could be issues with unique identifiers.

As with many of my blog posts, I’m writing this in the hope that it will appear in searches when someone else has the same problem that I ran into recently. In the past I’ve found myself Googling problems that I’ve had before and finding a link to my blog with an explanation in a post that I didn’t even remember writing.

Anyway, the problem I’m talking about today is one I ran into when working with a client a few weeks ago.

I was using the ogr2ogr command-line tool (part of the GDAL software suite) to import data from a local Geopackage file into a PostGIS database (ie. a PostgreSQL database with the PostGIS extension).

I had multiple files of data that I wanted to put into one Postgres table. Specifically, I was using the lovely data collated by Alasdair Rae on the resources page of his website. Even more specifically, I was using some of the Local Authority GIS data to get buildings data for various areas of the UK. I downloaded multiple GeoPackage files (for example, for Southampton City Council, Hampshire County Council and Portsmouth City Council) and wanted to import them all to a buildings table.

I originally tested this with a Postgres server running on my local machine, and ran the following ogr2ogr commands:

ogr2ogr --debug ON \
   -f PostgreSQL PG:"host=localhost user=postgres password=blah dbname=test_db" \
  buildings1.gpkg -nln buildings

ogr2ogr -append -update --debug ON \
   -f PostgreSQL PG:"host=localhost user=postgres password=blah dbname=test_db" \
   buildings2.gpkg -nln buildings

Here I’m using the -f switch and the arguments following it to tell ogr2ogr to export to PostgreSQL and how to connect to the server, giving it the input file of buildings1.gpkg and using the -nln parameter to tell it what layer name (ie. table name) to use as the output. In the second command I do exactly the same with buildings2.gpkg but also add -append and -update to tell it to append to the existing table rather than overwriting it.

This all worked fine. Great!

A few days later I tried the same thing with a Postgres server running on Azure (using Azure Database for PostgreSQL). The first command ran fine, but the second command seemed to hang.

I was expecting that it would be a bit slower when connecting to a remote database, but I left it running for 10 minutes and it still hadn’t finished. I then tried importing the second file to a new table and it completed quickly – therefore suggesting it was some sort of problem with appending the data.

I worked round this for the time being (using the ogrmerge.py script to merge my buildings1.gpkg and buildings2.gpkg into one file and then importing that file), but resolved to get to the bottom of it when I had time.

Recently, I had that time, and posted on the GDAL mailing list about this. The maintainer of GDAL got back to me to tell me about something I’d missed in the documentation. This was that when importing to a brand new table, the Postgres COPY mode is used, but when appending to an existing table individual INSERT statements are used instead, which can be a lot slower.

Let’s look into this in a bit more detail. The PostgreSQL COPY command is a fast way of importing data into Postgres which involves copying a whole file of data into Postgres in one go, rather than dealing with each row of data individually. This can be significantly faster than iterating through each row of the data and running a separate INSERT statement for each row.

So, ogr2ogr hadn’t hung, it was just running extremely slowly, as inserting my buildings layer involved running an INSERT statement separately for each row, and there were hundreds of thousands of rows. Because the server was hosted remotely on Azure, this involved sending the INSERT command from my computer to the server, waiting for the server to process it, and then the server sending back a result to my computer – a full round-trip for each row of the table.

So, I was told, the simple way to speed this up was to use a configuration setting to turn COPY mode on when appending to tables. This can be done by adding --config PG_USE_COPY YES to the ogr2ogr command. This did the job, and the append commands now completed nice and quickly. If you’re using GDAL/OGR from within a programming language, then have a look at the docs for the GDAL bindings for your language – there should be a way to set GDAL configuration options in your code.

The only final part of this was to understand why the COPY method isn’t used all the time, as it’s so much quicker. Even explained that this is because of potential issues with other connections to the database updating the table at the same time as GDAL is accessing it. It is a fairly safe assumption that if you’re creating a brand new table then no-one else will be accessing it yet, but you can’t assume the same for an existing table. The COPY mode can’t deal with making sure unique identifiers are unique when other connections may be accessing the data. whereas individual INSERT statements can cope with this. Therefore it’s safer to default to INSERT statements when there is any risk of data corruption.

As a nice follow-up for this, and on the maintainer’s advice, I submitted a PR to the GDAL docs, which adds a new section explaining this and giving guidance on setting the config option. I’ve copied that section below:

When data is appended to an existing table (for example, using the -append option in ogr2ogr) the driver will, by default, emit an INSERT statement for each row of data to be added. This may be significantly slower than the COPY-based approach taken when creating a new table, but ensures consistency of unique identifiers if multiple connections are accessing the table simultaneously.

If only one connection is accessing the table when data is appended, the COPY-based approach can be chosen by setting the config option PG_USE_COPY to YES, which may significantly speed up the operation.


How to create a Windows USB boot stick using Ventoy

I ran into a situation recently where I needed to create a Windows 10 bootable USB stick. I could easily download a Windows 10 ISO file, but I knew that it needed some ‘special magic’ to write it to a USB stick that would boot properly.

I tried various solutions (including windiskwriter) but none of them worked – until I tried Ventoy.

I got rather confused when reading the Ventoy webpage, so I wanted to write a quick blog post to show how to use it for this specific use case. The technology I had available to me was one M1 MacBook Pro running macOS and one x64 desktop machine with no OS (on which I wanted to install Windows) – and so these instructions will be based around this situation.

The way that Ventoy works is that you write the Ventoy image to a USB stick and it creates multiple partitions: some boot partitions, plus a big empty partition using up most of the space on the USB stick. You then copy your boot ISO files (for Windows, or Linux or whatever) onto that partition. When you boot from the USB stick Ventoy will provide a menu to allow you to pick an ISO, and then it will boot that ISO as if it had been burned to the USB stick itself. The easiest (and most reliable) way to write the Ventoy image to a USB stick is to use a bootable Ventoy image to do it – a double-USB-stick approach.

So, in my situation what I needed to do was:

  1. Download the Ventoy Live CD ISO file from https://ventoy.net/en/download.html
  2. Burn that to a USB stick (it can be a small one) using a standard USB stick burning program like balenaEtcher
  3. Boot from that USB stick on my x64 machine.
  4. Insert another USB stick (ideally a fairly large one) into that machine (while keeping the first one inserted) and follow the on-screen instructions to burn Ventoy to that USB stick.
  5. Shutdown that computer
  6. Insert that USB stick into my Mac, and copy the Windows 10 ISO file onto the large empty partition.
  7. Boot from that USB stick on the x64 computer. It will show an interface to allow you to choose an ISO to boot – there will be only one, so just press Enter. A few minutes later the Windows 10 installer will be ready to use.

Of course you can now keep this Ventoy USB stick around, and add other ISO files to it as needed.


How to fix GeoPandas drop_duplicates on geometry column not getting rid of all duplicates?

I use GeoPandas for a lot of my vector GIS data manipulation in Python.

I had a situation the other day where I ended up with duplicates of some geometries in my GeoDataFrame, and I wanted to remove them. The simple way to do this is to use the underlying pandas method drop_duplicates on the geometry column, like:

gdf.drop_duplicates('geometry')

However, even after running this I still had some duplicate geometries present. I couldn’t quite believe this, so checked multiple times and then started examining the geometries themselves.

What I found was that my duplicates were technically different geometries, but they looked the same when viewing them on a map. That’s because my geometries were LineStrings and I had two copies of the geometry: one with co-ordinates listed in the order left-to-right, and one in the order right-to-left.

This is illustrated in the image below: both lines look the same, but one line has the individual vertex co-ordinates in order from left-to-right and one has the same co-ordinates in order from right-to-left.

These two geometries will show as the same when using the geometry.equals() method, but won’t be picked up by drop_duplicates. That’s because drop_duplicates just serialises the geometry to Well-Known Binary and compares those to check for equality.

I started implementing various complex (and computationally-intensive) ways to deal with this, and then posted an issue on the GeoPandas Github page. Someone there gave me a simple solution which I want to share with you.

All you need to do is run gdf.normalize() first. So, the full code would be:

gdf.normalize()
gdf.drop_duplicates('geometry')

The normalize() method puts the vertices into a standard order so that they can be compared easily. This works for vertex order in lines and polygons, and ring orders in complex polygons.


How to easily reload a QGIS layer

Summary: The QGIS Reloader plugin will give you a new QGIS command to reload a layer, very useful when you’re updating a layer in another program and want to refresh QGIS’s view of it.

I use QGIS for almost all of my interactive GIS needs. Often I’ll be writing some code in Python that creates some sort of GIS layer (often a GeoPackage) and then I’ll use QGIS to interactively explore that layer.

In this situation, I’ll frequently create new versions of the GIS layer. Usually I’ve dealt with this by just creating XXX_v1.gpkg, XXX_v2.gpkg and so on, and dragging them into QGIS. However, if you have some complex symbology set up on the layer in QGIS you have the additional step of copying/pasting the layer styling from the old layer to the new layer. That’s a bit of a pain, but at least it’s possible.

Recently I’ve been using the QGIS Attribute Forms functionality to display an image when clicking on a feature in a layer, and those settings can’t be easily copied between layers. If I drag in a new version of the layer then I’ll have to click through and set up all of the form settings again.

What I wanted was a way to reload the contents of a layer from the file on disk – so I could just save to the same filename each time and reload the layer in QGIS. A bit of research showed me that duplicating a layer did this (right-click the layer and click Duplicate Layer) – but that leaves you with multiple layers, and an annoying ‘copy’ appended to the name.

Today I found a plugin that does exactly what I need. It’s called Reloader and is available from the QGIS Plugins Repository – just search for its name. It adds some new menu entries to reload the selected layers, reopen the layers (which updates the extent too), and even start/stop watching layers for changes (with a filesystem watcher that will update the layer whenever the underlying file changes). Just what I needed!

The only minor issue I have is that the menu entries come under the top-level Plugins menu, rather than in the right-click menu for an additional layer. I assume that QGIS doesn’t let plugins modify the layer context menu. However, to save going to the Plugins menu each time you can set up a keyboard shortcut for these functions. Go to Settings->Keyboard Shortcuts and search for Reload or Reopen and click the Change button to set a keyboard shortcut.

In the image above I’ve set a keyboard shortcut that looks very complex to press, but is actually easy when using my Goku configuration to get Caps Lock working as a modifier key.


How to create an x64 (Intel) conda environment on your Apple Silicon Mac (ARM) conda install

I came across some conda packages that didn’t work properly on my M1 Mac (Apple Silicon – ARM processor) the other day. They installed fine, but gave segmentation faults when run. So, I wanted to run the x64 (Intel) versions of these packages instead.

I haven’t actually needed to do this since I got a M1 Mac (a testament to the quality and range of Arm conda packages available these days through conda-forge), so I wasn’t sure how to do it.

A bit of Googling and experimenting led to this nice simple set of instructions. The upshot of this is you don’t need to install another version of Anaconda/miniconda/etc – you can just create a x64 environment in your existing conda install.

So:

  1. Run
    CONDA_SUBDIR=osx-64 conda create -n your_environment_name python

    This is a standard command to create a conda environment, but with CONDA_SUBDIR=osx-64 prepended to the command. This sets the CONDA_SUBDIR environment variable, and tells conda to use packages in the osx-64 subdirectory of the package server, rather than the standard osx-arm64 (M1 Mac) subdirectory. This is what gets you the x64 packages.

  2. When this command finishes, you will have a new x64 conda environment. You can then activate it with
    conda activate your_environment_name
  3. Now we need to tell conda to always use this CONDA_SUBDIR setting when using this environment, otherwise any future installs in this environment will use the default CONDA_SUBDIR and things will get very confused. We can do this by setting a conda environment config setting to tell conda to set a specific environment variable when you activate the environment. Do this by running:
    conda env config vars set CONDA_SUBDIR=osx-64
  4. The output of that command will warn you to deactivate and reactivate the environment, so do this
    conda deactivate
    conda activate your_environment_name
  5. That’s it! You now have a x64 environment which you can install packages in using standard conda install commands.