Category Archives: qgis

Properly setting up your QGIS license

If you want your copy of QGIS display it’s legal licensing status, this is the missing code for you.

Copy this in your QGIS Python script editor (WARNING: DO NOT RUN THIS IN AN IMPORTANT USER PROFILE, I will NOT help you if it breaks something):

import os
from qgis.core import QgsApplication, QgsSettings
from qgis.PyQt.QtGui import QColor, QImage, QPainter, QPen, QStaticText
from qgis.PyQt.QtWidgets import QMessageBox

def install_qgis_license():
    # WARNING: This fucks around with your profiles and stuff!
    # QgsCustomization is not available from Python so just yolo here and write a fresh file if possible
    profile_directory = QgsApplication.qgisSettingsDirPath()
    customization_ini_filepath = profile_directory + "QGIS/QGISCUSTOMIZATION3.ini"

    if os.path.isfile(customization_ini_filepath):
        # ain't gonna be touchin dat!
        text = (
            "QGISCUSTOMIZATION3.ini EXISTS! UNLICENSED HACKING DETECTED!\n"
            "Or: Custom license has been installed already...\n"
            "Anyways, we are not creating a *new* license now ;)"
        )
        messagebox = QMessageBox()
        messagebox.setText(text);
        messagebox.exec()
        return

    # get existing splash image
    splash_path = QgsApplication.splashPath()  # :/images/splash/
    splash_image_file = splash_path + "splash.png"
    splash_image = QImage(splash_image_file)

    # paint new splash image
    new_splash_image = QImage(splash_image)
    painter = QPainter()
    painter.begin(new_splash_image)

    # white bold font plz
    font = QgsApplication.font()
    font.setBold(True)
    painter.setFont(font)
    pen = painter.pen()
    pen.setColor(QColor("white"))
    painter.setPen(pen)

    # place text at appropriate location
    label_text = f"This QGIS©®™ is legally licensed to {os.getlogin()}"
    label_text_rect = painter.boundingRect(0, 0, 0, 0, 0, label_text)  # returns new rect that fits text
    left_offset = new_splash_image.width() - label_text_rect.size().width() - 20
    painter.drawText(left_offset, 840, label_text)

    painter.end()

    # create license dir if necessary
    new_splash_dir_path = profile_directory + "license"
    try:
        os.mkdir(new_splash_dir_path)
    except FileExistsError:
        pass

    save_success = new_splash_image.save(new_splash_dir_path + "/splash.png")
    if save_success:
        print(f"Initialized new QGIS license....")
    else:
        print("Error on QGIS license initialization, this will get reported!")
        return

    # enable license dir for splash image lookup in QGISCUSTOMIZATION3.ini
    with open(customization_ini_filepath, "w") as sink:
        sink.write("[Customization]\n")
        sink.write(f"splashpath={new_splash_dir_path}/")  # must have trailing slash

    # enable loading of QGISCUSTOMIZATION3.ini in user profile
    QgsSettings().setValue(r"UI/Customization/enabled", True)
    
    print("License installed, reboot QGIS to activate!")
    messagebox = QMessageBox()
    messagebox.setText("License installed, restart QGIS now to activate!");
    messagebox.exec()

#install_qgis_license()

Then (if you really want to do it), uncomment the function call in the last line and execute the script. Follow the instructions.

To clean up remove or restore the QGIS/QGISCUSTOMIZATION3.ini file in your profile and remove the license directory from your profile, restore the previous value of UI/Customization/enabled in your profile (just remove the line or disable Settings -> Interface Customization).

If you want to hate yourself in the future, put it in a file called startup.py in QStandardPaths.standardLocations(QStandardPaths.AppDataLocation) aka the directory which contains the profiles directory itself.

BTW: If you end up with QGIS crashing and lines like these in the error output:

...
Warning: QPaintDevice: Cannot destroy paint device that is being painted
QGIS died on signal 11
...

It is probably not a Qt issue that caused the crash. The QPaintDevice warning might just be Qt telling you about your painter being an issue during clean up of the actual crash (which might just be a wrong name or indentation somewhere in your code, cough).

Animated snowflakes in QGIS

Inspired by https://mapstodon.space/@bogind/109558968869019746 and listening to https://www.youtube.com/watch?v=AdyAF9M3XVw on repeat thanks to https://ruby.social/@halfbyte/109558803178815394 I made:

This was the final expression (with lots of opportunity to improve):

with_variable(
  'point_at_top_of_canvas',
  densify_by_count(
    make_line(
      point_n( @map_extent, 3),  -- no idea if these indexes are stable
      point_n( @map_extent, 4)
    ),
    42  -- number of trajectories
  ),
  collect_geometries(
    array_foreach(
      generate_series(1, num_points(@point_at_top_of_canvas)),
      with_variable(
        'point_n_of_top_line',
        point_n(@point_at_top_of_canvas, @element),
        point_n(
          wave_randomized(
            make_line(
              @point_n_of_top_line,
              -- make it at least touch the bottom of the canvas:
              translate(@point_n_of_top_line, 0, -@map_extent_height)
            ),
            -- fairly stupid frequency and wavelength but hey, works in any crs
            1, @map_extent_width/5,
            1, @map_extent_width/100,
            seed:=@element  -- stable waves \o/
          ),
          floor(epoch(now())%10000/50)  -- TODO make it loop around according to num_points of each line
        )
      )
    )
  )
)

Use it on an empty polygon layer with an inverted polygon style and set it to refresh at a high interval (0.01s?). Or use this QGIS project (I included some intermediate steps of the style as layer styles if you want to learn about this kind of stuff):

#30DayMapChallenge as #1Day30MapsChallenge (2021)

Not sure why I never posted this last year but I did the #30DayMapChallenge in a single day, streamed live via a self-hosted Owncast instance. It was … insane and fun. This year I will do it again, on the 26th of November.

Here are most of the maps I made last year:

Some notes I kept, please bug me about recovering the others from my Twitter archive (I deleted old tweets a bit too early):

  • 1 Points: Pins via Geometry Generator in QGIS
  • 2 Lines: River elevation profiles of Elbe, Rhein, Ems, Weser, Donau and Main. DEM: © GeoBasis-DE / BKG (2021)
  • 13 NaEr I mean Natural Earth (Blame @tjukanov)
  • 18 Water (DGM-W 2010 Unter- und Außenelbe, Wasserstraßen- und Schifffahrtsverwaltung des Bundes, http://kuestendaten.de, 2010)
  • 20 Movement: Emojitions on a curvy trajectory. State changes depending on the curvyness ahead. Background: (C) OpenStreetMap Contributors <3
  • 21 Elevation with qgis2threejs (It’s art, I swear!
  • 22 Boundaries: Inspired by Command and Conquer Red Alert. Background by Spiney (CC-BY 3.0 / CC-BY-SA 3.0, https://opengameart.org/node/12098)
  • 24 Historical: Buildings in Hamburg that were built before the war (at least to some not so great dataset). Data Lizenz: Datenlizenz Deutschland Namensnennung 2.0 (Freie und Hansestadt Hamburg, Landesbetrieb Geoinformation und Vermessung (LGV))
  • 27 Heatmap: Outdoor advertisements (or something like that) in Hamburg. Fuck everything about that! Data Lizenz: Datenlizenz Deutschland Namensnennung 2.0 (Freie und Hansestadt Hamburg, Behörde für Verkehr und Mobilitätswende, (BVM))
  • 28 Earth not flat. Using my colleague’s Beeline plugin to create lines between the airports I have flown too and the Globe Builder plugin by @gispofinland to make a globe.

Waggawaggawaggawagga animated ducks in QGIS

I used ne_110m_admin_0_countries.

Rendering updates for the layer at 0.1 seconds.

Geometry Generator for a Point for the marker location via line_interpolate_point:

with_variable(
  'biggest_geom',
  geometry_n(order_parts($geometry, 'area($geometry)', ascending:=False), 1),
  line_interpolate_point(
    boundary(@biggest_geom), 
    perimeter(@biggest_geom)*(round(epoch(now())/100)%100/100)
  )
)

Raster Image Marker with https://opengameart.org/content/character-spritesheet-duck, vertical anchor at bottom, sprite choice between walking and running (doesn’t actually work) plus the frame via

with_variable(
  'biggest_geom',
  geometry_n(order_parts($geometry, 'area($geometry)', ascending:=False), 1),
  '/your/path/Duck/Sprites/Walking-Running/'
  || if(perimeter(@biggest_geom) < 10, 'Walking', 'Running')
  || ' 00'
  || to_string(round(epoch(now())/200)%2+1)
  || '.png'
)

Rotation did not work, I tried line_interpolate_angle:

with_variable(
  'biggest_geom',
  geometry_n(order_parts($geometry, 'area($geometry)', ascending:=False), 1),
  line_interpolate_angle(
    boundary(@biggest_geom), 
    perimeter(@biggest_geom)*(round(epoch(now())/100)%100/100)
  )
)

Steps via two more Geometry Generators, both for Lines using line_substring and some nice style (inspired by the wonderful built-in cat trail preset):

with_variable(
	'biggest_geom',
	geometry_n(order_parts($geometry, 'area($geometry)', ascending:=False), 1),
  	line_substring(
	  boundary(@biggest_geom), 
	  0,
	  perimeter(@biggest_geom)*(round(epoch(now())/100)%100/100)
	)
)

Could be improved if (for example) Raster Image Marker would support:

  • Choice of resampling algorithm
  • Flipping
  • Rotation would work, no idea what’s wrong with my expression, it works with random values, so …
  • Whatever is broken with the choice between ‘Walking’ and ‘Running’ in the file path expression

Highlight current timeslice in a QGIS Atlas layout

Did this for an ex-colleague some months ago and forgot to share the how-to publically. We needed a visual representation of the current time in a layout that showed both a raster map (different layer per timeslice) and a timeseries plot of an aspect of the data (this was created outside QGIS).

Have lots of raster layers you want to iterate through. I have:

./ECMWF_ERA_40_subset/2019-01-01.tif
./ECMWF_ERA_40_subset/2019-01-02.tif
./ECMWF_ERA_40_subset/2019-01-03.tif
./ECMWF_ERA_40_subset/2019-01-04.tif
./ECMWF_ERA_40_subset/2019-01-05.tif
./ECMWF_ERA_40_subset/2019-01-06.tif
...

Create a new layer for your map extent. Draw your extent as geometry. Duplicate that geometry as many times as you have days. Alternatively you could of course have different geometries per day. Whatever you do, you need a layer with one feature per timeslice for the Atlas to iterate though. I have 30 days to visualise so I duplicated my extent 30 times.

Open the Field Calculator. Add a new field called date as string type (not as date type until some bug is fixed (sorry, did not make a note here, maybe sorting is/was broken?)) with an expression that represents time and orders chronologically if sorted by QGIS. For example: '2019-01-' || lpad(@row_number,2,0) (assuming your records are in the correct order if you have different geometries…)

Have your raster layers named the same way as the date attribute values.

Make a new layout.

For your Layout map check “Lock layers” and use date as expression for the “Lock layers” override. This will now select the appropriate raster layer, based on the attribute value <-> layer name, to display for each Atlas page.

Cool, if you preview the Atlas now you got a nice animation through your raster layers. Let’s do part 2:

In your layout add your timeseries graph. Give it a unique ID, e. g. “plot box”. Set its width and height via new variables (until you can get those via an expression this is needed for calculations below).

Create a box to visualise the timeslice. Set its width to map_get(item_variables('plot box'), 'plot_width') / @atlas_totalfeatures. For the height and y use/adjust this expression: map_get(item_variables('plot box'), 'plot_height'). For x comes the magic:

with_variable(
'days_total',
day(to_date(maximum("date"))-to_date(minimum("date")))+1,
-- number of days in timespan
-- +1 because we need the number of days in total
-- not the inbetween, day() to just get the number of days
with_variable(
'mm_per_day',
map_get(item_variables('plot box'), 'plot_width') / @days_total,
with_variable(
'days',
day(to_date(attribute(@atlas_feature, 'date'))-to_date(minimum("date"))),
-- number of days the current feature is from the first day
-- to_date because BUG attribute() returns datetime for date field
@mm_per_day * @days + map_get(item_variables('plot box'), 'plot_x')
)
)
)

This will move the box along the x axis accordingly.

Have fun!

How to get a dataset out of OpenStreetMap (OSM) and into QGIS easily

For small datasets just use the QuickOSM plugin, enter your key=value and it will load the data directly as QGIS layer(s).

For big datasets don’t run massive queries on the Overpass API but use prepared thematic or regional OpenStreetMap extracts.

For example, if you want all things tagged place=town globally, grab place_EPSG4326.zip from http://download.osmdata.xyz/ and run a filter on that locally.

Don’t be scared of processing a global OSM planet file either, Osmium makes filtering OSM data extremely easy:

osmium tags-filter in.osm.pbf n/place=town -o place=town.osm.pbf

Yes, a planet file is pretty big, but extracting specific features from that is not a big data problem and you must not be scared of it. Downloading the file will probably take you magnitudes longer than extracting something from it. For me it was 45 minutes for the download, then about 8 minutes for extracting on a seriously slow (~60MB/s) spinning metal hard disk drive (no SSD).

And yes, you can load OSM PBF directly into QGIS thanks to GDAL’s support for the format.

While we are at it: Don’t use GeoJSON for anything but data transfer and maybe storage. It is not an efficient format to power your layers and GIS analyses. (OSM PDF isn’t either.)

#30DayMapChallenge: Day 3 – Polygons OR Lego® style brick raster in QGIS using Geometry Generator expressions

I showed this at the annual meeting of the swiss user group in June 2019 and promptly forgot to post it for everyone to see. Let’s blame the scenery?

Have a raster.

Processing -> Create Grid, covering it with points in a spacing of your choice. Use the same CRS as your raster (unless you want to figure out expression-based geometry transformations on your own, like in my elevation lines code)

Change the Symbol layer type to Geometry Generator and enter

with_variable(
  'radius',
  3333,
  buffer($geometry, @radius, 16)
)

where radius should be a value about one third of your spacing.

You should see circles!

For the fill color use this expression and adjust the name of your raster layer:

with_variable(
    'raster_layer',
    'DHM200.xyz',
    ramp_color('Blues',  -- change to an other named ramp here if you like
        scale_linear(
            raster_value(
                @raster_layer, 1, 
                centroid($geometry)  -- back to our point
            ),
            raster_statistic(
                @raster_layer, 1, 
                'min'
            ),
            raster_statistic(
                @raster_layer, 1, 
                'max'
            ),
            0, 1  -- new scale as color ramps go 0 to 1
        )
    )
)

This will get the raster value below our grid point and fit it onto a color ramp between the min and max of all the raster values.

For the stroke use the same expression but wrap darker(..., 150) around it so you get a darker color.

Using Draw Effects add a small Drop Shadow to your circles. I used an offset of 0.5 mm and a blur radius of 1 mm.

Now add another Geometry Generator symbol layer below your existing one and use the following expression:

bounds(buffer($geometry, 5000))

with_variable(
  'radius',
  5000,
  bounds(
    buffer($geometry, @radius)
  )
)

with the radius being half your grid spacing.

Use the same expressions for the colors as above but set the darker value to 200.

For some more fanciness maybe add a “QGIS” text on top of the nupsies?

Exercise: Make it so that the result perfectly covers the raster, instead of being one grid cell off like mine.

Thanks Topi.

#30DayMapChallenge: Day 2 – Lines

Thanks Topi.

Addresses in Hamburg: http://suche.transparenz.hamburg.de/dataset/alkis-adressen-hamburg15

Convert the housenumber from string to integer: to_int("hausnr")

Processing -> Points to Path using above housenumber as order field and strname as group field.

Color those ziggyzaggety lines as you like.

Maybe explode the lines and style each segment by its housenumber?

And put them on an interactive slippy map.