Advertisements

Archive

QGIS

I’ll start with some tech talk first. Feel free to jump to the usage example further down if you are here for the edge bundling plugin.

As you certainly know, QGIS 3 brings a lot of improvements and under-the-hood changes. One of those changes affects all Python scripts. They need to be updated to Python 3 and the new PyQGIS API. (See the official migration guide for details.)

To get ready for the big 3.0 release, I’ve started porting my Processing tools. The edge bundling script is my first candidate for porting to QGIS 3. I also wanted to use this opportunity to “upgrade” from a simple script to a plugin that integrates into Processing.

I used Alexander Bruy’s “prepair for Processing” plugin as a template but you can also find an example template in your Processing folder. (On my system, it is located in C:\OSGeo4W64\apps\qgis-dev\python\plugins\processing\algs\exampleprovider.)

Since I didn’t want to miss the advantages of a good IDE, I set up PyCharm as described by Heikki Vesanto. This will give you code completion for Python 3 and PyQGIS which is very helpful for refactoring and porting. (I also tried Eclipse with PyDev but if you don’t have a favorite IDE yet, I find PyCharm easier to install and configure.)

My PyCharm startup script qgis3_pycharm.bat is a copy of C:\OSGeo4W64\bin\python-qgis-dev.bat with the last line altered to start PyCharm:

@echo off
call "%~dp0\o4w_env.bat"
call qt5_env.bat
call py3_env.bat
@echo off<span data-mce-type="bookmark" style="display: inline-block; width: 0px; overflow: hidden; line-height: 0;" class="mce_SELRES_start"></span>
path %OSGEO4W_ROOT%\apps\qgis-dev\bin;%PATH%
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT:\=/%/apps/qgis-dev
set GDAL_FILENAME_IS_UTF8=YES
rem Set VSI cache to be used as buffer, see #6448
set VSI_CACHE=TRUE
set VSI_CACHE_SIZE=1000000
set QT_PLUGIN_PATH=%OSGEO4W_ROOT%\apps\qgis-dev\qtplugins;%OSGEO4W_ROOT%\apps\qt5\plugins
set PYTHONPATH=%OSGEO4W_ROOT%\apps\qgis-dev\python;%PYTHONPATH%
start /d "C:\Program Files\JetBrains\PyCharm\bin\" pycharm64.exe

In PyCharm File | Settings, I configured the OSGeo4W Python 3.6 interpreter and added qgis-dev and the plugin folder to its path:

With this setup done, we can go back to the code.

I first resolved all occurrences of import * in my script to follow good coding practices. For example:

from qgis.core import *

became

from qgis.core import QgsFeature, QgsPoint, QgsVector, QgsGeometry, QgsField, QGis<span data-mce-type="bookmark" style="display: inline-block; width: 0px; overflow: hidden; line-height: 0;" class="mce_SELRES_start"></span>

in this PR.

I didn’t even run the 2to3 script that is provided to make porting from Python 2 to Python 3 easier. Since the edge bundling code is mostly Numpy, there were almost no changes necessary. The only head scratching moment was when Numpy refused to add a map() return value to an array. So (with the help of Stackoverflow of course) I added a work around to convert the map() return value to an array as well:

flocal_x = map(forcecalcx, subtr_x, subtr_y, distance)
electrostaticforces_x[e_idx, :] += np.array(list(flocal_x))

The biggest change related to Processing is that the VectorWriter has been replaced by a QgsFeatureSink. It’s defined as a parameter of the edgebundling QgsProcessingAlgorithm:

self.addParameter(QgsProcessingParameterFeatureSink(
   self.OUTPUT,
   self.tr("Bundled edges"),
   QgsProcessing.TypeVectorLine)
)

And when the algorithm is run, the sink is filled with the output features:

(sink, dest_id) = self.parameterAsSink(
   parameters, self.OUTPUT, context,
   source.fields(), source.wkbType(), source.sourceCrs()
)

# code that creates features

sink.addFeature(feat, QgsFeatureSink.FastInsert)

The ported plugin is available on Github.

The edge bundling plugin in action

I haven’t uploaded the plugin to the official plugin repository yet, but you can already download if from Github and give it a try:

For this example, I’m using taxi pick-up and drop-off data provided by the NYC Taxi & Limousine Commission. I downloaded the January 2017 green taxi data and extracted all trips for the 1st of January. Then I created origin-destination (OD) lines using the QGIS virtual layer feature:

To get an interesting subset of the data, I extracted only those OD flows that cross the East River and have a count of at least 5 taxis:

Now the data is ready for bundling.

If you have installed the edge bundling plugin, the force-directed edge bundling algorithm should be available in the Processing toolbox. The UI of the edge bundling algorithm looks pretty much the same as it did for the QGIS 2 Processing script:

Since this is a small dataset with only 148 OD flows, the edge bundling processes is pretty quick and we can explore the results:

Beyond this core edge bundling algorithm, the repository also contains two more scripts that still need to be ported. They include dependencies on sklearn, so it will be interesting to see how straightforward it is to convert them.

Advertisements

QGIS 3 has a new feature: reports! In short, reports are the good old Altas feature on steroids.

Let’s have a look at an example project:

To start a report, go to Project | New report. The report window is quite similar to what we’ve come to expect from Print Composer (now called Layouts). The most striking difference is the report panel at the left side of the screen.

When a new report is created, the center of the report window is empty. To get started, we need to select the report entry in the panel on the left. By selecting the report entry, we get access to the Include report header and Include report footer checkboxes. For example, pressing the Edit button next to the Include report header option makes it possible to design the front page (or pages) of the report:

Similarly, pressing Edit next to the Include report footer option enables us to design the final pages of our report.

Now for the content! We can populate our report with content by clicking on the plus button to add a report section or a “field group”. A field group is basically an Atlas. For example, here I’ve added a field group that creates one page for each country in the Natural Earth countries layer that I have loaded in my project:

Note that in the right panel you can see that the Controlled by report option is activated for the map item. (This is equivalent to a basic Atlas setup in QGIS 2.x.)

With this setup, we are ready to export our report. Report | Export Report as PDF creates a 257 page document:

As configured, the pages are ordered by country name. This way, for example, Australia ends up on page 17.

Of course, it’s possible to add more details to the individual pages. In this example, I’ve added an overview map in Robinson projection (to illustrate again that it is now possible to mix different CRS on a map).

Happy QGIS mapping!

If you have already designed a few maps in QGIS, you are probably aware of a long-standing limitation: Print Composer maps were limited to the project’s coordinate reference system (CRS). It was not possible to have maps with different CRS in a composition.

Note how I’ve been using the past tense? 

Rejoice! QGIS 3 gets rid of this limitation. Print Composer has been replaced by the new Layout dialog which – while very similar at first sight – offers numerous improvements. But today, we’ll focus on projection handling.

For example, this is a simple project using WGS84 as its project CRS:


In the Layouts dialog, each map item now has a CRS property. For example, the overview map is set to World_Robinson while the main map is set to ETRS-LAEA:

As you can see, the red overview frame in the upper left corner is curved to correctly represent the extent of the main map.

Of course, CRS control is not limited to maps. We also have full freedom to add map grids in yet another CRS:

This opens up a whole new level of map design possibilities.

Bonus fact: Another great improvement related to projections in QGIS3 is that Processing tools are now aware of layers with different CRS and will actively reproject layers. This makes it possible, for example, to intersect two layers with different CRS without any intermediate manual reprojection steps.

Happy QGIS mapping!

Data Plotly is a new plugin by Matteo Ghetta for QGIS3 which makes it possible to draw D3 graphs of vector layer attribute values. This is a huge step towards making QGIS a one stop shop for data exploration!

Data Plotly adds a new panel where graphs can be configured and viewed. Currently, there are nine different plot types:

The following examples use tree cadastre data from the city of Linz, Austria.

Scatter plots with both two and three variables are supported. After picking the attributes you want to visualize, press “Create plot”.

If you change some settings and press “Create plot” again, by default, the new graph will be plotted on top of the old one. If you don’t want that to happen, press “Clean plot canvas” before creating a new plot.

The plots are interactive and display more information on mouse over, for example, the values of a box plot:

Even aggregate expressions are supported! Here’s the mean height of trees by type (deciduous L or coniferous N):

For more examples, I strongly recommend to have a look at the plugin home page.

In this post, I want to show how to visualize building block data published by the city of Vienna in 3D using QGIS. This data is interesting due to its level of detail. For example, here you can see the Albertina landmark in the center of Vienna:

an this is the corresponding 3D visualization, including flying roof:

To enable 3D view in QGIS 2.99 (soon to be released as QGIS 3), go to View | New 3D Map View.

Viennese building data (https://www.data.gv.at/katalog/dataset/76c2e577-268f-4a93-bccd-7d5b43b14efd) is provided as Shapefiles. (Saber Razmjooei recently published a similar post using data from New York City in ESRI Multipatch format.) You can download a copy of the Shapefile and a DEM for the same area from my dropbox.  The Shapefile contains the following relevant attributes for 3D visualization

  • O_KOTE: absolute building height measured to the roof gutter(?) (“absolute Gebäudehöhe der Dachtraufe”)
  • U_KOTE: absolute height of the lower edge of the building block if floating above ground (“absolute Überbauungshöhe unten”)
  • HOEHE_DGM: absolute height of the terrain (“absolute Geländehöhe”)
  • T_KOTE: lowest point of the terrain for the given building block (“tiefster Punkt des Geländes auf den Kanten der Gebäudeteilfläche”)

To style the 3D view in QGIS 3, I set height to “U_KOTE” and extrusion to

O_KOTE-coalesce(U_KOTE,0)

both with a default value of 0 which is used if the field or expression is NULL:

The altitude clamping setting defines how height values are interpreted. Absolute clamping is perfect for the Viennese data since all height values are provided as absolute measures from 0. Other options are “relative” and “terrain” which add given elevation values to the underlying terrain elevation. According to the source of qgs3dutils:

  AltClampAbsolute,   //!< Z_final = z_geometry
  AltClampRelative,   //!< Z_final = z_terrain + z_geometry
  AltClampTerrain,    //!< Z_final = z_terrain

The gray colored polygon style shown in the map view on the top creates the illusion of shadows in the 3D view:

 

Beyond that, this example also features elevation model data which can be configured in the 3D View panel. I found it helpful to increase the terrain tile resolution (for example to 256 px) in order to get more detailed terrain renderings:

Overall, the results look pretty good. There are just a few small glitches in the rendering, as well as in the data. For example, the kiosik in front of Albertina which you can also see in the StreetView image, is lacking height information and therefore we can only see it’s “shadow” in the 3D rendering.

So far, I found 3D rendering performance very good. It works great on my PC with Nvidia graphics card. On my notebook with Intel Iris graphics, I’m unfortunately still experiencing crashes which I hope will be resolved in the future.

MarineCadastre.gov is a great source for AIS data along the US coast. Their data formats and tools though are less open. Luckily, GDAL – and therefore QGIS – can read ESRI File Geodatabases (.gdb).

MarineCadastre.gov also offer a Track Builder script that creates lines out of the broadcast points. (It can also join additional information from the vessel and voyage layers.) We could reproduce the line creation step using tools such as Processing’s Point to path but this post will show how to create PostGIS trajectories instead.

First, we have to import the points into PostGIS using either DB Manager or Processing’s Import into PostGIS tool:

Then we can create the trajectories. I’ve opted to create a materialized view:

The first part of the query creates a temporary table called ptm (short for PointM). This step adds time stamp information to each point. The second part of the query then aggregates these PointMs into trajectories of type LineStringM.

CREATE MATERIALIZED VIEW ais.trajectory AS
 WITH ptm AS (
   SELECT b.mmsi,
     st_makepointm(
       st_x(b.geom), 
       st_y(b.geom), 
       date_part('epoch', b.basedatetime)
     ) AS pt,
     b.basedatetime t
   FROM ais.broadcast b
   ORDER BY mmsi, basedatetime
 )
 SELECT row_number() OVER () AS id,
   st_makeline(ptm.pt) AS st_makeline,
   ptm.mmsi,
   min(ptm.t) AS min_t,
   max(ptm.t) AS max_t
 FROM ptm
 GROUP BY ptm.mmsi
WITH DATA;

The trajectory start and end times (min_t and max_t) are optional but they can help speed up future queries.

One of the advantages of creating trajectory lines is that they render many times faster than the original points.

Of course, we end up with some artifacts at the border of the dataset extent. (Files are split by UTM zone.) Trajectories connect the last known position before the vessel left the observed area with the position of reentry. This results, for example, in vertical lines which you can see in the bottom left corner of the above screenshot.

With the trajectories ready, we can go ahead and start exploring the dataset. For example, we can visualize trajectory speed and/or create animations:

Purple trajectory segments are slow while green segments are faster

We can also perform trajectory analysis, such as trajectory generalization:

This is a first proof of concept. It would be great to have a script that automatically fetches the datasets for a specified time frame and list of UTM zones and loads them into PostGIS for further processing. In addition, it would be great to also make use of the information in the vessel and voyage tables, thus splitting up trajectories into individual voyages.


This post is part of a series. Read more about movement data in GIS.

Today’s post is a follow-up of Movement data in GIS #3: visualizing massive trajectory datasets. In that post, I summarized a concept for trajectory generalization. Now, I have published the scripts and sample data in my QGIS-Processing-tools repository on Github.

To add the trajectory generalization scripts to your Processing toolbox, you can use the Add scripts from files tool:

It is worth noting, that Add scripts from files fails to correctly import potential help files for the scripts but that’s not an issue this time around, since I haven’t gotten around to actually write help files yet.

The scripts are used in the following order:

  1. Extract characteristic trajectory points
  2. Group points in space
  3. Compute flows between cells from trajectories

The sample project contains input data, as well as output layers of the individual tools. The only required input is a layer of trajectories, where trajectories have to be LINESTRINGM (note the M!) features:

Trajectory sample based on data provided by the GeoLife project

In Extract characteristic trajectory points, distance parameters are specified in meters, stop duration in seconds, and angles in degrees. The characteristic points contain start and end locations, as well as turns and stop locations:

The characteristic points are then clustered. In this tool, the distance has to be specified in layer units, which are degrees in case of the sample data.

Finally, we can compute flows between cells defined by these clusters:

Flow lines scaled by flow strength and cell centers scaled by counts

If you use these tools on your own data, I’d be happy so see what you come up with!


This post is part of a series. Read more about movement data in GIS.

%d bloggers like this: