Archive

Author Archives: underdark

Update 2017: Please follow the instructions on the OSGeo4W website: https://trac.osgeo.org/osgeo4w/wiki/ExternalPythonPackages


Update 2015: I now recommend this workflow based on Nathan’s post on installing Python setup tools:

  1. Download ez_setup.py
  2. Run python ez_setup.py in your OSGeo4W shell
  3. Run easy_install pysal in your OSGeo4W shell

Original post:

Today’s post is a summary of how to install PySAL on Windows for OSGeo4W Python.

The most important resource was https://trac.osgeo.org/osgeo4w/wiki/ExternalPythonPackages

In the OSGeo4W Shell run:

C:\Users\anita_000\Desktop>curl http://python-distribute.org/distribute_setup.py | python

Update: Note that python-distribute.org has gone down since I posted this. See http://www.gossamer-threads.com/lists/python/python/1158958 for more info.

Then download https://raw.githubusercontent.com/pypa/pip/master/contrib/get-pip.py to the Desktop and run:

C:\Users\anita_000\Desktop>python get-pip.py

When pip is ready, install PySAL:

C:\Users\anita_000\Desktop>pip install pysal

Test the installation:

C:\Users\anita_000\Desktop>python
Python 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)] on win 32
Type "help", "copyright", "credits" or "license" for more information.
>>> import pysal

The point table of the Spatialite database created from OSM north-eastern Austria contains more than 500,000 points. This post shows how the style works which – when applied to the point layer – wil make sure that only towns and (when zoomed in) villages will be marked and labeled.

Screenshot 2014-07-12 12.30.21

In the attribute table, we can see that there are two tags which provide context for populated places: the place and the population tag. The place tag has it’s own column created by ogr2ogr when converting from OSM to Spatialite. The population tag on the other hand is listed in the other_tags column.

Screenshot 2014-07-12 13.00.15

for example

"opengeodb:lat"=>"47.5000237","opengeodb:lon"=>"16.0334769","population"=>"623"

Overview maps would be much too crowded if we simply labeled all cities and towns. Therefore, it is necessary to filter towns based on their population and only label the bigger ones. I used limits of 5,000 and 10,000 inhabitants depending on the scale.

Screenshot 2014-07-12 12.56.33

At the core of these rules is an expression which extracts the population value from the other_tags attribute: The strpos() function is used to locate the text "population"=>" within the string attribute value. The population value is then extracted using the left() function to get the characters between "population"=>" and the next occurrence of ". This value can ten be cast to integer using toint() and then compared to the population limit:

5000 < toint( 
   left (
      substr(
         "other_tags",
         strpos("other_tags" ,'"population"=>"')+16,
         8
      ),
      strpos(
         substr(
            "other_tags",
            strpos("other_tags" ,'"population"=>"')+16,
            8
         ),
        '"'
      )
   )
) 

There is also one additional detail concerning label placement in this style: When zoomed in closer than 1:400,000 the labels are placed on top of the points but when zoomed out further, the labels are put right of the point symbol. This is controlled using a scale-based expression in the label placement:

Screenshot 2014-07-12 13.32.47

As usual, you can find the style on Github: https://github.com/anitagraser/QGIS-resources/blob/master/qgis2/osm_spatialite/osm_spatialite_tonerlite_point.qml

In my opinion, Stamen’s Toner-lite map is one of the best background maps to use together with colorful overlays. The only downsides of using it in QGIS are that the OpenLayers plugin can not provide the tiles at print resolution and that the projection is limited to Web Mercator. That’s why I’ve started to recreate the style for OSM Spatialite databases:

toner-lite

So far, there are styles for lines and polygons and they work quite well for the scale range between 1:1 and 1:250000. As always, you can download the styles from QGIS-resources on Github.

Using OSM data in QGIS is a hot topic but so far, no best practices for downloading, preprocessing and styling the data have been established. There are many potential solutions with all their advantages and disadvantages. To give you a place to start, I thought I’d share a workflow which works for me to create maps like the following one from nothing but OSM:

osm_google_100k

Getting the data

Raw OSM files can be quite huge. That’s why it’s definitely preferable to download the compressed binary .pbf format instead of the XML .osm format.

As a download source, I’d recommend Geofabrik. The area in the example used in this post is part of the region Pays de la Loire, France.

Preparing the data for QGIS

In the preprocessing step, we will extract our area of interest and convert the .pbf into a spatialite database which can be used directly in QGIS.

This can be done in one step using ogr2ogr:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat 2.59 46.58 -1.44 47.07 noirmoutier.db noirmoutier.pbf

where the -spat option controls the area of interest to be extracted.

When I first published this post, I suggested a two step approach. You can find it here for future reference:

For the first step: extracting the area of interest, we need Osmosis. (For Windows, you can get osmosis from openstreetmap.org. Unpack to use. Requires Java.)

When you have Osmosis ready, we can extract the area of interest to the .osm format:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>..\bin\osmosis.bat --read-pbf pays-de-la-loire-latest.osm.pbf --bounding-box left=-2.59 bottom=46.58 right=-1.44 top=47.07 --write-xml noirmoutier.osm

While QGIS can also load .osm files, I found that performance and access to attributes is much improved if the .osm file is converted to spatialite. Luckily, that’s easy using ogr2ogr:

C:\Users\anita_000\Geodata\OSM_Noirmoutier>ogr2ogr -f "SQLite" -dsco SPATIALITE=YES noirmoutier.db noirmoutier.osm

Finishing preprocessing in QGIS

In QGIS, we’ll want to load the points, lines, and multipolygons using Add SpatiaLite Layer:

Screenshot 2014-05-31 11.39.40

When we load the spatialite tables, there are a lot of features and some issues:

  • There is no land polygon. Instead, there are “coastline” line features.
  • Most river polygons are missing. Instead there are “riverbank” line features.

Screenshot 2014-05-31 11.59.58

Luckily, creating the missing river polygons is not a big deal:

  1. First, we need to select all the lines where waterway=riverbank.
    Screenshot 2014-05-31 13.14.00
  2. Then, we can use the Polygonize tool from the processing toolbox to automatically create polygons from the areas enclosed by the selected riverbank lines. (Note that Processing by default operates only on the selected features but this setting can be changed in the Processing settings.)
    Screenshot 2014-05-31 13.40.16

Creating the land polygon (or sea polygon if you prefer that for some reason) is a little more involved since most of the time the coastline will not be closed for the simple reason that we are often cutting a piece of land out of the main continent. Therefore, before we can use the Polygonize tools, we have to close the area. To do that, I suggest to first select the coastline using "other_tags" LIKE '%"natural"=>"coastline"%' and create a new layer from this selection (save selection as …) and edit it (don’t forget to enable snapping!) to add lines to close the area. Then polygonize.

Screenshot 2014-05-31 14.38.48

Styling the data

Now that all preprocessing is done, we can focus on the styling.

You can get the styles used in the map from my Github QGIS-resources repository:

  • osm_spatialite_googlemaps_multipolygon.qml … rule-based renderer incl. rules for: water, natural, residential areas and airports
  • osm_spatialite_googlemaps_lines.qml … rule-based renderer incl. rules for roads, rails, and rivers, as well as rules for labels
  • osm_spatialite_googlemaps_roadshields.qml … special label style for road shields
  • osm_spatialite_googlemaps_places.qml … label style for populated places such as cities and towns

qgis_osm_google_100k

This post shows how to quickly and easily create a small QGIS plugin for counting the number of features within a vector layer.

To get started, you will need QGIS and Qt Designer (to design the user interface) installed. If you are on Windows, I suggest WinPython which provides Qt Designer and Spyder (a Python IDE).

The great thing about creating plugins for QGIS: There is a plugin for that! It’s called Plugin Builder. And while you are at it, also install Plugin Reloader. Reloader is very useful for plugin developers because it lets you quickly reload your plugin without having to restart QGIS every time you make changes to the code.

installPluginBuilder

Plugin Builder will create all the files we need for our plugin. Just start it and select a name for your plugin class (one word in CamelCase), as well as a name for the plugin itself and the plugin menu entry (can be multiple words). Once you press Ok, you’re asked to select a folder to store the plugin. You can save directly to the QGIS plugin folder ~\.qgis2\python\plugins.

pluginBuilder

Next, open the newly created folder (in my case ~\.qgis2\python\plugins\BuilderTest). Amongst other files, it contains the user interface file ui_buildertest.ui. Our plugin will count the number of features in a vector layer. Therefore, it needs a combobox which allows the user to select a layer. Open the .ui file in Qt Designer and add a combobox to the dialog. Change the object name of the combobox to layerCombo. We’ll later use this name in the plugin code to add items to the combobox. Save the dialog and close Qt Designer.

qtDesigner

Now, we need to compile the .ui and the resources.qrc file to turn the dialog and the icon into usable Python code. This is done on the command line. On Windows, I suggest using the OSGeo4W Shell. Navigate to the plugin folder and run:

pyuic4 -o ui_buildertest.py ui_buildertest.ui
pyrcc4 -o resources_rc.py resources.qrc

Note: Using the latest version of Plugin Builder, you only need to compile resources.qrc because the .ui is now loaded dynamically. Furthermore, you should use the following command to compile the resources:

pyrcc4 -o resources.py resources.qrc 

If you enable and run the plugin now, you will already see the dialog but the combobox will be empty. To populate the combobox, we need to write a few lines of code in buildertest.py. First, we’ll fetch all loaded layers and add all vector layers to the combobox. Then, we’ll add code to compute and display the number of features in the selected layer. To achieve this, we expand the run() method:

def run(self):        
    # show the dialog
    self.dlg.show()

    layers = QgsMapLayerRegistry.instance().mapLayers().values()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.layerCombo.addItem( layer.name(), layer ) 
         
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result == 1:
        # do something useful 
        index = self.dlg.layerCombo.currentIndex()
        layer = self.dlg.layerCombo.itemData(index)
        QMessageBox.information(self.iface.mainWindow(),"hello world","%s has %d features." %(layer.name(),layer.featureCount()))

When you are done with the code, you can use Plugin Reloader to load the new version. When you start the plugin now, the combobox will be populated with the names of the vector layers in your current project. And on pressing Ok, the plugin will compute and display the number of features.

builderTEst

builderTestResult

For more information on PyQGIS and more code samples I warmly recommend the PyQGIS Cookbook. Have fun!

Extracting POIs from OpenStreetMap is reasonably simple using Overpass API. A very convenient way to construct the query is to use a query builder which allows you to select the area of interest and builds queries for different servers.

xapi_query_builder

Of course you can fine-tune the query further. For example, you can add multiple key-value pairs to the query. I used the following query to select all Billa supermarkets:

www.overpass-api.de/api/xapi?*[shop=supermarket][name=Billa][bbox=15.96725,48.0432,16.79947,48.40915]

Note the * in the query? It means that I’m querying all kinds of features: nodes, ways, and relations.

Save the server response to a .osm file. This file can be loaded into QGIS using simple drag-and-drop or Add Vector Layer. A dialog will open where you can select the type of features you want to load from the file. You can simply use Select All and OK to load everything.

add_osm_file

My supermarket POIs came in two types: points and multipolygons. To style them both with nice supermarket SVG icons, I decided to use a Centroid fill with the SVG marker for the polygon layer:

osm_pois_billa

Open data and open source GIS … nice :-)

Today was the last day of the Vienna code sprint which brought together OSGeo developers from many projects. It’s been a great week thanks to organizers and sponsors!

The QGIS team was extremely busy working on the project’s web infrastructure (e.g. new plugins.qgis.org website) as well as hunting down and fixing bugs.

Check out some impressions on twitter.

qgis23_vienna

More pictures on the official blog: vienna2014.sprint.osgeo.org

If you are looking for a tool to easily create 3D visualizations of your geodata, look no further! Qgis2threejs is a plugin by Minoru Akagi which exports terrain data combined with the map canvas image and optional vector data to an html file which can be viewed in 3D in any web browser which supports WebGL. To do that, this plugin uses the Three.js library.

This is the result of my first experiments with Qgis2threejs. In the following sections, I will show the steps to reproduce it.

Türkenschanzpark, Vienna

click for the interactive version (requires WebGL-capable browser)

1. The data

The building blocks of this visualization are:

  • elevation data and the hillshade derived from this data
  • a base map (WMTS from basemap.at in my case)
  • OSM building data provided by Geofabrik and
  • tree data from the city of Vienna

Load all datasets into QGIS.

2. Preparing the map

Qgis2threejs will overlay the map (as rendered in the QGIS map area) on top of the elevation model. You can combine any number of layers to create your map. I just loaded a basemap.at WMTS and a hillshade layer. To add a nice tree shadow effect, I also added the tree layer (dark grey, 50% transparency, multiply blending).

tuerkenschanzpark_map

3. Preparing the vector features

The vector features in the visualization are buildings and trees. The buildings are based on an OSM building layer. The trees are create from two point layers: one point layer to create the tree trunks (cylinder shape) and a duplicate of this point layer to create the tree crowns (sphere shape).

Load the data and choose the desired fill colors.

4. Using Qgis2threejs

Now we can start Qgis2threejs. The first tab is used to configure the terrain. Just pick the correct elevation data layer. I didn’t modify any of the other default settings.

qgis2threejs_dem

The second tab provides the settings for the vector data. As mentioned in the previous section, the trees are created from two point layers and the buildings are based on a polygon layer. The tree crowns are spheres with a radius size 3 and a z value of 5 above the surface. The tree trunks are cylinders. Finally, the buildings have a height of 10.

qgis2threejs_vector

That’s it! Just press “run” and wait. When the export is finished, your default browser (or a different one, if you specify another one in the plugin settings) will open automatically and display the results.
The visualization is interactive. You can tilt the visualization using the left mouse button, pan using the right mouse button, and zoom using the mouse wheel. I found that Firefox used around 1.6 GB of RAM to render this example.

5. Share your visualization

In the browser window, you will see where Qgis2threejs stored the html and associated Javascript files. To share your visualization, you just need to copy these files onto a webserver.

I would love to see what you come up with. Please share a link in the comments.

QGIS 2.2 is now available for Windows through OSGeo4W installer. Packages for other systems are being prepared by the package maintainers.

The Windows packages are currently marked experimental, so you have to use the advanced install in OSGeo4W and check the ‘Exp’ radio button on the top to install them.

osgeo4w_qgis22

As release manager Jürgen Fischer announced:

Please test and report problems, so that I can soon promote them to ‘curr’ent.
Once that has happend, I’ll proceed with turning them into standalone
installers.

QGIS 2.2 will be released tomorrow, February 21st. Following the release of 2.0, the QGIS project decided to move to a time-based release plan with releases every four months. This provides a clear framework for developers, translators and documenters which makes it possible to plan ahead and know when tasks have to be finished to be included in a release version.

Similar to the 2.0 release, the project has invested considerable resources to make 2.2 “Valmiera” a successful release. I have already blogged about some of the great new features. Thanks to the project donors and sponsors it was also possible to fund developer time for many important bug fixes.

One of the greatest resources of the QGIS project are its users. One group that deserves our special thanks is the Swiss QGIS User Group. They collect a modest annual membership fee which provides a steady and growing crowd-funding that can be used to positively influence the QGIS project. For example, they invested in bug fixing for 2.0 and they are co-funding work on multi-threaded rendering for QGIS 2.4.

With the rise of new QGIS user groups “QUGs” all around the world, e.g. in Australia, the UK, and the US, I hope these groups will find ways to bring users together and to positively influence the development of QGIS towards the next releases.