Related Plugins and Tags

QGIS Planet

GRASS GIS 6.4.3RC4 released

Fourth (and last) release candidate of GRASS GIS 6.4.3 with improvements and stability fixes
A fourth release candidate of GRASS GIS 6.4.3 is now available.

Source code download:

Binaries download:

To get the GRASS GIS 6.4.3RC4 source code directly from SVN:
 svn checkout http://svn.osgeo.org/grass/grass/tags/release_20130710_grass_6_4_3RC4

Key improvements of this release include some new functionality (assistance for topologically unclean vector data), fixes in the vector network modules, fixes for the wxPython based portable graphical interface (attribute table management, wxNVIZ, and Cartographic Composer), fixes in the location wizard for Datum transform selection and support for PROJ.4 version 4.8.0, improvements for selecting the Python version to be used, enhanced portability for MS-Windows (native support, fixes in case of missing system DLLs), and more translations (esp. Romanian).

See also our detailed announcement:
 http://trac.osgeo.org/grass/wiki/Release/6.4.3RC4-News

First time users should explore the first steps tutorial after installation.

Release candidate management at
http://trac.osgeo.org/grass/wiki/Grass6Planning

Please join us in testing this release candidate for the final release.

Consider to donate pizza or beer for the upcoming GRASS GIS Community Sprint in Prague:
Thanks to all contributors!

Support the upcoming GRASS GIS Community Sprint with a (micro)donation!

Have you been using GRASS GIS within your company or for your daily tasks? Would like to express your satisfaction and give something back to these terrific software developers to make them even happier during the July’s GRASS GIS Community Sprint in Prague?
Sure – so don’t think twice and support the developers with a (micro)donation!
Please contact Markus Neteler ([email protected] – GRASS GIS PSC Chair) for further details.
Wait – it is also possible to buy a round of beer for the developers with a quick click using the PayPal “Buy [pizza/beer/...] Now” button below: (… more pizza & beer for the developers)

About the Community Sprint

A “Community Sprint” is a get-together for members and supporters of GRASS GIS and related OSGeo projects to make decisions and tackle larger problems. Developers and contributors are donating their valuable time, so appreciate direct or in-kind funding made available for the sprint meeting to cover out-of-pocket expenses. All of the work that takes place at the community sprint will be directly contributed back into the GRASS GIS project to the benefit of everyone who uses it.
See the outstanding results from 2011, and 2012!

Public transport isochrones with pgRouting

This post covers a simple approach to calculating isochrones in a public transport network using pgRouting and QGIS.

For this example, I’m using the public transport network of Vienna which is loaded into a pgRouting-enable database as network.publictransport. To create the routable network run:

select pgr_createTopology('network.publictransport', 0.0005, 'geom', 'id');

Note that the tolerance parameter 0.0005 (units are degrees) controls how far link start and end points can be apart and still be considered as the same topological network node.

To create a view with the network nodes run:

create or replace view network.publictransport_nodes as
select id, st_centroid(st_collect(pt)) as geom
from (
	(select source as id, st_startpoint(geom) as pt
	from network.publictransport
	) 
union
	(select target as id, st_endpoint(geom) as pt
	from network.publictransport
	) 
) as foo
group by id;

To calculate isochrones, we need a cost attribute for our network links. To calculate travel times for each link, I used speed averages: 15 km/h for buses and trams and 32km/h for metro lines (similar to data published by the city of Vienna).

alter table network.publictransport add column length_m integer;
update network.publictransport set length_m = st_length(st_transform(geom,31287));

alter table network.publictransport add column traveltime_min double precision;
update network.publictransport set traveltime_min = length_m  / 15000.0 * 60; -- average is 15 km/h
update network.publictransport set traveltime_min = length_m  / 32000.0 * 60 where "LTYP" = '4'; -- average metro is 32 km/h

That’s all the preparations we need. Next, we can already calculate our isochrone data using pgr_drivingdistance, e.g. for network node #1:

create or replace view network.temp as
 SELECT seq, id1 AS node, id2 AS edge, cost, geom
  FROM pgr_drivingdistance(
    'SELECT id, source, target, traveltime_min as cost FROM network.publictransport',
    1, 100000, false, false
  ) as di
  JOIN network.publictransport_nodes pt
  ON di.id1 = pt.id;

The resulting view contains all network nodes which are reachable within 100,000 cost units (which are minutes in our case).

Let’s load the view into QGIS to visualize the isochrones:

isochrone_publictransport_1

The trick is to use data-defined size to calculate the different walking circles around the public transport stops. For example, we can set up 10 minute isochrones which take into account how much time was used to travel by pubic transport and show how far we can get by walking in the time that is left:

1. We want to scale the circle radius to reflect the remaining time left to walk. Therefore, enable Scale diameter in Advanced | Size scale field:

scale_diameter

2. In the Simple marker properties change size units to Map units.
3. Go to data defined properties to set up the dynamic circle size.

datadefined_size

The expression makes sure that only nodes reachable within 10 minutes are displayed. Then it calculates the remaining time (10-"cost") and assumes that we can walk 100 meters per minute which is left. It additionally multiplies by 2 since we are scaling the diameter instead of the radius.

To calculate isochrones for different start nodes, we simply update the definition of the view network.temp.

While this approach certainly has it’s limitations, it’s a good place to start learning how to create isochrones. A better solution should take into account that it takes time to change between different lines. While preparing the network, more care should to be taken to ensure that possible exchange nodes are modeled correctly. Some network links might only be usable in one direction. Not to mention that there are time tables which could be accounted for ;)


pgRouting 2.0 for Windows quick guide

This post is a quick instruction for installing Postgres 9.2, PostGIS 2.0 and pgRouting 2.0.

  1. For Postgres, download the installer from enterprisedb.com.
  2. Run the installer. You’ll have to pick a superuser password – remember it, you’ll need it again soon.
  3. At the end of the installation process, allow the installer to start Stack Builder.
  4. In Stack Builder, select the Postgres 9.2 installation and install PostGIS from the list of available extensions.
  5. The PostGIS installation procedure will prompt you for the superuser password you picked before.
  6. I suggest letting the installer create a sample database We’ll need it later anyway.

Now for the pgRouting part:

  1. Download the pgRouting zip file for your system (32 or 64 bit) from Winnie.
  2. Unzip the file. It contains bin, lib and share folders as well as two text files.
  3. Copy these folders and files over to your Postgres installation. In my case C:\Program Files\PostgreSQL\9.2

Installation – done.

Next, fire up pgAdmin. If you allowed the PostGIS installer to create a sample database, you will find it named postgis20 or similar. Otherwise just create a new database using the PostGIS template database. You can enable pgRouting in a database using the following steps:

  1. In postgis20, execute the following query to create the pgrouting extension. This will add the pgRouting-specific functions:
    CREATE EXTENSION pgrouting;
  2. Test if everything worked fine:
    SELECT pgr_version();

    It should return "(2.0.0-dev,v2.0.0-beta,18,a3be38b,develop,1.46.1)" or similar – depending on the version you downloaded.

pgadmin_pgrouting

How about some test data? I’ll be using the public transport network of the city of Vienna provided in GeoJSON format from http://data.wien.gv.at/daten/wfs?service=WFS&request=GetFeature&version=1.1.0&typeName=ogdwien:OEFFLINIENOGD&srsName=EPSG:4326&outputFormat=json:

  1. Just copy paste the url in Add Vector Layer | Protocol to load the dataset.
  2. Use DB Manager to load the layer into your database. (As you can see in the screenshot, I created a schema called network but that’s optional.)
  3. import_publictransport

  4. To make the line vector table routable, we use pgr_createTopology. This function assumes the columsn called “source” and “target” exist so we have to create those as well:
    alter table network.publictransport add column source integer;
    alter table network.publictransport add column target integer;
    select pgr_createTopology('network.publictransport', 0.0001, 'geom', 'id');
    

    I’m quite generously using a tolerance of 0.0001 degrees to build the topology. Depending on your dataset, you might want to be more strict here.

  5. Let’s test it! Route from source #1 to target #3000 using pgr_dijkstra:
    SELECT seq, id1 AS node, id2 AS edge, cost, geom
      FROM pgr_dijkstra(
        'SELECT id, source, target, st_length(geom) as cost FROM network.publictransport',
        1, 3000, false, false
      ) as di
      JOIN network.publictransport pt
      ON di.id2 = pt.id ;

    Note how the query joins the routing results and the network table together. (I’m aware that using the link length as a cost attribute will not lead to realistic results in a public transport network but bear with me for this example.)

  6. We can immediately see the routing results using the Load as new layer option:

select_dijkstra
route

Nice work! pgRouting 2.0 has come a long way. In a post from April this year, Boston GIS even announced to add pgRouting into the Stack Builder. That’s going to make the installation even more smooth.


OpenLayers for QGIS master

As mentioned in my previous post, all plugins need an update to function in QGIS master – soon to be 2.0. The OpenLayers plugin is one of the most popular plugins but so far has not been updated by the developer. If you miss it, you can get a fixed version from qgis.nl Temporary Fix for OpenLayers Plugin thanks to Richard!

openlayers


Temporary Fix for OpenLayers Plugin

Current development version of QGIS ( upcoming QVIS 2.0 ) has undergone a nessecary upgrade of the Python-Cpp glue (SIP) which temporarily broke almost all python plugins. Currently plugin devs are busy fixing their plugins to be usable on time for the real launch of QGIS 2.0. One plugin notably missing plugin is the OpenLayers-plugin, […]

Alpha by Value choropleth in QGIS

This post was inspired by a post written by Xaquín G.V, some extra credit must also go to Mathieu for helping me with the expression.

One of the new great features in the upcoming QGIS 2.0 is the ability to use data defined symbology and labelling. It might look like nothing on the surface but it can really open up the possibility of some great looking maps.

Lets start with the example Xaquín has used by classifying by the unemployment rate. I have even created a pre-joined dataset for you if you want to follow along, grab it here.

So first lets make a custom colour ramp just for this:

usaramp.png

call it "usa" and classify the data

usaclass.png

usaresult1.png

Ok not bad. But you will always have those hard edges on the colours. What if we could get the exact colour ramp value for a given value? Well you can.

We need to make a new column for the colour value for each feature. Open up the field calculator and add a new column with the following expression

ramp_color('usa', scale_linear( "unemployed_by_county_xgv_Rate",0,18,0,1 ))

What are the ramp_colur and scale_linear functions? ramp_colur is a new function that allows you to get a RGBA value for a position inside the ramp. So if you want the RGBA value at the 75% mark you can use ramp_color('usa', 0.75).

Cool but to make it dynamic we need to bring in the value from the "Rate" column, however that goes up higher then 1. How do we get in a range from 0..1? Well we can scale it back down using scale_linear.

scale_linear(<column>,<min in value>,<max in value>,<min out value>,<max out value>)

Having said all that here is our result:

usaattribute.png

To use this new colour column we need to use data defined symbols

usasymbol.png

Note: I have selected Single Symbol

And the result

usaresult2.png

This isn't going to suit every map but when something is scale based this method works well.

Bumping it up a notch

Higher unemployment rates plus high lobor force really should get more weight. One way to handle this is by using Alpha-By-Value like Xaquín said in his post. Well to do that in QGIS we need to change the alpha value of what is returned from color_ramp. Time for a bit of regex (don't worry in 2.1 I will add a colourramprgb function with no aplha part)

First lets get just the aplha value part by scaling the lobor force coloum. Add a new column using the following expression

toint(scale_linear("unemployed_by_county_xgv_Labor_Force", 0, 100000, 0, 255))

usaattributealhpa.png

Now make a new column that replaces the alpha part of the colour column

regexp_replace( "colour", ',[^,]*$' , format(',%1',"alpha" ))

usaattributealphacolour.png

Then we can use data defined symbols on that new column

usaresult3.png

Nifty!

Now you will have to play with the min and max values for the scaling functions and also the colours in the ramp but you get the idea.

Bonus

You can also use the same colour column on the labels to match the polygon and even increase the label size based on the rate column

usalabels.png

So that is your Thursday lesson on how to create a Alpha by value map in the upcoming QGIS 2.0.

Extra stuff

If you want to avoid creating all those extra columns or even do it on the fly you can use this single expression

regexp_replace( ramp_color('usa', scale_linear( "unemployed_by_county_xgv_Rate",0,15,0,1)),',[^,]*$',','|| toint(scale_linear("unemployed_by_county_xgv_Labor_Force",0,100000,0,255)))

Using the single expression block can be handy if you want to get your upper and lower limits right.

Alpha by Value choropleth in QGIS

This post was inspired by a post written by Xaquín G.V, some extra credit must also go to Mathieu for helping me with the expression.

One of the new great features in the upcoming QGIS 2.0 is the ability to use data defined symbology and labelling. It might look like nothing on the surface but it can really open up the possibility of some great looking maps.

Lets start with the example Xaquín has used by classifying by the unemployment rate. I have even created a pre-joined dataset for you if you want to follow along, grab it here.

So first lets make a custom colour ramp just for this:

Alt Text

call it "usa" and classify the data

Alt Text

Alt Text

Ok not bad. But you will always have those hard edges on the colours. What if we could get the exact colour ramp value for a given value? Well you can.

We need to make a new column for the colour value for each feature. Open up the field calculator and add a new column with the following expression

ramp_color('usa', scale_linear( "unemployed_by_county_xgv_Rate",0,18,0,1 ))

What are the ramp_colur and scale_linear functions? ramp_colur is a new function that allows you to get a RGBA value for a position inside the ramp. So if you want the RGBA value at the 75% mark you can use ramp_color('usa', 0.75).

Cool but to make it dynamic we need to bring in the value from the "Rate" column, however that goes up higher then 1. How do we get in a range from 0..1? Well we can scale it back down using scale_linear.

scale_linear(<column>,<min in value>,<max in value>,<min out value>,<max out value>) 

Having said all that here is our result:

Alt Text

To use this new colour column we need to use data defined symbols

Alt Text

Note: I have selected Single Symbol

And the result

Alt Text

This isn't going to suit every map but when something is scale based this method works well.

Bumping it up a notch

Higher unemployment rates plus high lobor force really should get more weight. One way to handle this is by using Alpha-By-Value like Xaquín said in his post. Well to do that in QGIS we need to change the alpha value of what is returned from color_ramp. Time for a bit of regex (don't worry in 2.1 I will add a colour_ramp_rgb function with no aplha part)

First lets get just the aplha value part by scaling the lobor force coloum. Add a new column using the following expression

toint(scale_linear("unemployed_by_county_xgv_Labor_Force", 0, 100000, 0, 255))

Alt Text

Now make a new column that replaces the alpha part of the colour column

regexp_replace( "colour", ',[^,]*$' , format(',%1',"alpha" ))

Alt Text

Then we can use data defined symbols on that new column

Alt Text

Nifty!

Now you will have to play with the min and max values for the scaling functions and also the colours in the ramp but you get the idea.

Bonus

You can also use the same colour column on the labels to match the polygon and even increase the label size based on the rate column

Alt Text

So that is your Thursday lesson on how to create a Alpha by value map in the upcoming QGIS 2.0.

Extra stuff

If you want to avoid creating all those extra columns or even do it on the fly you can use this single expression

regexp_replace( ramp_color('usa', scale_linear( "unemployed_by_county_xgv_Rate",0,15,0,1)),',[^,]*$',','|| toint(scale_linear("unemployed_by_county_xgv_Labor_Force",0,100000,0,255)))

Using the single expression block can be handy if you want to get your upper and lower limits right.

Alpha by Value choropleth in QGIS

This post was inspired by a post written by Xaquín G.V, some extra credit must also go to Mathieu for helping me with the expression.

One of the new great features in the upcoming QGIS 2.0 is the ability to use data defined symbology and labelling. It might look like nothing on the surface but it can really open up the possibility of some great looking maps.

Lets start with the example Xaquín has used by classifying by the unemployment rate. I have even created a pre-joined dataset for you if you want to follow along, grab it here.

So first lets make a custom colour ramp just for this:

Alt Text

call it "usa" and classify the data

Alt Text

Alt Text

Ok not bad. But you will always have those hard edges on the colours. What if we could get the exact colour ramp value for a given value? Well you can.

We need to make a new column for the colour value for each feature. Open up the field calculator and add a new column with the following expression

ramp_color('usa', scale_linear( "unemployed_by_county_xgv_Rate",0,18,0,1 ))

What are the ramp_colur and scale_linear functions? ramp_colur is a new function that allows you to get a RGBA value for a position inside the ramp. So if you want the RGBA value at the 75% mark you can use ramp_color('usa', 0.75).

Cool but to make it dynamic we need to bring in the value from the "Rate" column, however that goes up higher then 1. How do we get in a range from 0..1? Well we can scale it back down using scale_linear.

scale_linear(<column>,<min in value>,<max in value>,<min out value>,<max out value>) 

Having said all that here is our result:

Alt Text

To use this new colour column we need to use data defined symbols

Alt Text

Note: I have selected Single Symbol

And the result

Alt Text

This isn't going to suit every map but when something is scale based this method works well.

Bumping it up a notch

Higher unemployment rates plus high lobor force really should get more weight. One way to handle this is by using Alpha-By-Value like Xaquín said in his post. Well to do that in QGIS we need to change the alpha value of what is returned from color_ramp. Time for a bit of regex (don't worry in 2.1 I will add a colour_ramp_rgb function with no aplha part)

First lets get just the aplha value part by scaling the lobor force coloum. Add a new column using the following expression

toint(scale_linear("unemployed_by_county_xgv_Labor_Force", 0, 100000, 0, 255))

Alt Text

Now make a new column that replaces the alpha part of the colour column

regexp_replace( "colour", ',[^,]*$' , format(',%1',"alpha" ))

Alt Text

Then we can use data defined symbols on that new column

Alt Text

Nifty!

Now you will have to play with the min and max values for the scaling functions and also the colours in the ramp but you get the idea.

Bonus

You can also use the same colour column on the labels to match the polygon and even increase the label size based on the rate column

Alt Text

So that is your Thursday lesson on how to create a Alpha by value map in the upcoming QGIS 2.0.

Extra stuff

If you want to avoid creating all those extra columns or even do it on the fly you can use this single expression

regexp_replace( ramp_color('usa', scale_linear( "unemployed_by_county_xgv_Rate",0,15,0,1)),',[^,]*$',','|| toint(scale_linear("unemployed_by_county_xgv_Labor_Force",0,100000,0,255)))

Using the single expression block can be handy if you want to get your upper and lower limits right.

Adding __geo_interface__ to QGIS geometry and feature

What is __geo_interface__ anyway? If you have never seen __geo_interface__ it was defined by Sean Gillies at https://gist.github.com/sgillies/2217756. A cool lightweight interface that describes a single spatial object using a GeoJSON like structure, in fact parts of it are just normal GeoJSON.

Since QGIS in feature freeze at the moment I can't add this to the main code base just yet. Awwww sad face :(

but wait!

Thanks to good ol' Python we can just monkey patch it right on.

def mapping_feature(feature):
    geom = feature.geometry()
    properties = {}
    fields = [field.name() for field in feature.fields()]
    properties = dict(zip(fields, feature.attributes()))
    return { 'type' : 'Feature',
             'properties' : properties,
             'geometry' : geom.__geo_interface__}

def mapping_geometry(geometry):
    geo = geometry.exportToGeoJSON()
    # We have to use eval because exportToGeoJSON() gives us
    # back a string that looks like a dictionary. 
    return eval(geo)

QgsFeature.__geo_interface__ = property(lambda self: mapping_feature(self))
QgsGeometry.__geo_interface__ = property(lambda self: mapping_geometry(self))

And that's that. Easy hey. Really got to love dynamic languages.

The __geo_interface__ property now exists on any instance of QgsGeometry or QgsFeature. Lets test that theory.

>>> import pprint
>>> layer = iface.activeLayer()
>>> feature = layer.selectedFeatures()[0]
>>> feature.__geo_interface__
>>> pprint.pprint(f.__geo_interface__)
{'geometry': {'coordinates': [[[385039.90567724, 6449154.61878853],
                               [385059.01135993, 6449154.80874077],
                               [385059.41538719, 6449114.58680192],
                               [385040.30145863, 6449114.38685169],
                               [385040.16953133, 6449127.46423076],
                               [385039.90567724, 6449154.61878853]]],
              'type': 'Polygon'},
 'properties': {u'address': u'HYNES WY',
                u'assessment': u'2204315',
                u'bool': u'F',
                u'dola_pin': 283678,
                u'field18': 2920,
                u'field19': 0.0,
                u'house_numb': u'3',
                u'location': None,
                u'lot': u'LOT 107',
                u'new': None,
                u'old': u'http://www.seabreeze.com.au|mylink',
                u'paid_in_fu': u'F',
                u'pin_string': u'283678',
                u'postcode': u'6163',
                u'reserve': None,
                u'sample_dat': u'2003-05-15',
                u'subdivided': None,
                u'suburb': u'HAMILTON HILL',
                u'suburb_wit': u'HAMILTON HILL',
                u'type': u'H',
                u'v_auto_dat': None,
                u'v_auto_use': None,
                u'v_boolean': None,
                u'v_datetime': None,
                u'v_decimal': None,
                u'v_int': None,
                u'v_numeric': None,
                u'v_varcha_1': None,
                u'v_varchar': None,
                u'v_varchar_': None},
 'type': 'Feature'}
>>> feature.geometry().__geo_interface__
{'type': 'Point', 'coordinates': [388197.74503284, 6450504.16670842]}

Excellent!

Adding __geo_interface__ to QGIS geometry and feature

What is __geo_interface__ anyway? If you have never seen __geo_interface__ it was defined by Sean Gillies at https://gist.github.com/sgillies/2217756. A cool lightweight interface that describes a single spatial object using a GeoJSON like structure, in fact parts of it are just normal GeoJSON.

Since QGIS in feature freeze at the moment I can't add this to the main code base just yet. Awwww sad face :(

but wait!

Thanks to good ol' Python we can just monkey patch it right on.

def mapping_feature(feature):
    geom = feature.geometry()
    properties = {}
    fields = [field.name() for field in feature.fields()]
    properties = dict(zip(fields, feature.attributes()))
    return { 'type' : 'Feature',
             'properties' : properties,
             'geometry' : geom.__geo_interface__}

def mapping_geometry(geometry):
    geo = geometry.exportToGeoJSON()
    # We have to use eval because exportToGeoJSON() gives us
    # back a string that looks like a dictionary. 
    return eval(geo)

QgsFeature.__geo_interface__ = property(lambda self: mapping_feature(self))
QgsGeometry.__geo_interface__ = property(lambda self: mapping_geometry(self))

And that's that. Easy hey. Really got to love dynamic languages.

The __geo_interface__ property now exists on any instance of QgsGeometry or QgsFeature. Lets test that theory.

>>> import pprint
>>> layer = iface.activeLayer()
>>> feature = layer.selectedFeatures()[0]
>>> feature.__geo_interface__
>>> pprint.pprint(f.__geo_interface__)
{'geometry': {'coordinates': [[[385039.90567724, 6449154.61878853],
                               [385059.01135993, 6449154.80874077],
                               [385059.41538719, 6449114.58680192],
                               [385040.30145863, 6449114.38685169],
                               [385040.16953133, 6449127.46423076],
                               [385039.90567724, 6449154.61878853]]],
              'type': 'Polygon'},
 'properties': {u'address': u'HYNES WY',
                u'assessment': u'2204315',
                u'bool': u'F',
                u'dola_pin': 283678,
                u'field18': 2920,
                u'field19': 0.0,
                u'house_numb': u'3',
                u'location': None,
                u'lot': u'LOT 107',
                u'new': None,
                u'old': u'http://www.seabreeze.com.au|mylink',
                u'paid_in_fu': u'F',
                u'pin_string': u'283678',
                u'postcode': u'6163',
                u'reserve': None,
                u'sample_dat': u'2003-05-15',
                u'subdivided': None,
                u'suburb': u'HAMILTON HILL',
                u'suburb_wit': u'HAMILTON HILL',
                u'type': u'H',
                u'v_auto_dat': None,
                u'v_auto_use': None,
                u'v_boolean': None,
                u'v_datetime': None,
                u'v_decimal': None,
                u'v_int': None,
                u'v_numeric': None,
                u'v_varcha_1': None,
                u'v_varchar': None,
                u'v_varchar_': None},
 'type': 'Feature'}
>>> feature.geometry().__geo_interface__
{'type': 'Point', 'coordinates': [388197.74503284, 6450504.16670842]}

Excellent!

Adding __geo_interface__ to QGIS geometry and feature

What is __geo_interface__ anyway? If you have never seen __geo_interface__ it was defined by Sean Gillies at https://gist.github.com/sgillies/2217756. A cool lightweight interface that describes a single spatial object using a GeoJSON like structure, in fact parts of it are just normal GeoJSON.

Since QGIS in feature freeze at the moment I can't add this to the main code base just yet. Awwww sad face :(

but wait!

Thanks to good ol' Python we can just monkey patch it right on.

def mapping_feature(feature):
    geom = feature.geometry()
    properties = {}
    fields = [field.name() for field in feature.fields()]
    properties = dict(zip(fields, feature.attributes()))
    return { 'type' : 'Feature',
             'properties' : properties,
             'geometry' : geom.__geo_interface__}

def mapping_geometry(geometry):
    geo = geometry.exportToGeoJSON()
    # We have to use eval because exportToGeoJSON() gives us
    # back a string that looks like a dictionary. 
    return eval(geo)

QgsFeature.__geo_interface__ = property(lambda self: mapping_feature(self))
QgsGeometry.__geo_interface__ = property(lambda self: mapping_geometry(self))

And that's that. Easy hey. Really got to love dynamic languages.

The __geo_interface__ property now exists on any instance of QgsGeometry or QgsFeature. Lets test that theory.

>>> import pprint
>>> layer = iface.activeLayer()
>>> feature = layer.selectedFeatures()[0]
>>> feature.__geo_interface__
>>> pprint.pprint(f.__geo_interface__)
{'geometry': {'coordinates': [[[385039.90567724, 6449154.61878853],
                               [385059.01135993, 6449154.80874077],
                               [385059.41538719, 6449114.58680192],
                               [385040.30145863, 6449114.38685169],
                               [385040.16953133, 6449127.46423076],
                               [385039.90567724, 6449154.61878853]]],
              'type': 'Polygon'},
 'properties': {u'address': u'HYNES WY',
                u'assessment': u'2204315',
                u'bool': u'F',
                u'dola_pin': 283678,
                u'field18': 2920,
                u'field19': 0.0,
                u'house_numb': u'3',
                u'location': None,
                u'lot': u'LOT 107',
                u'new': None,
                u'old': u'http://www.seabreeze.com.au|mylink',
                u'paid_in_fu': u'F',
                u'pin_string': u'283678',
                u'postcode': u'6163',
                u'reserve': None,
                u'sample_dat': u'2003-05-15',
                u'subdivided': None,
                u'suburb': u'HAMILTON HILL',
                u'suburb_wit': u'HAMILTON HILL',
                u'type': u'H',
                u'v_auto_dat': None,
                u'v_auto_use': None,
                u'v_boolean': None,
                u'v_datetime': None,
                u'v_decimal': None,
                u'v_int': None,
                u'v_numeric': None,
                u'v_varcha_1': None,
                u'v_varchar': None,
                u'v_varchar_': None},
 'type': 'Feature'}
>>> feature.geometry().__geo_interface__
{'type': 'Point', 'coordinates': [388197.74503284, 6450504.16670842]}

Excellent!

Combining the best of Top10NL, BAG & OSM in QGIS

For some time now I’ve been using the BAG buildings in QGIS as an alternative to the Top10NL buildings. The BAG buildings are more accurate, after all. This is not to disqualify the Top10NL, on the contrary: each geo data set has its own merits. The Top10NL is unbeatable regarding terrain details. As for the […]

Raster Data Extraction using QIS

Raster files consist of a grid of cells, each cell contains a numeric value, which is used to determine how to colour each cell.  This value may be based on the elevation of the cell, flood water depth, or soil quality. It is possible to extract this information by point sampling or using a terrain profile. Point sampling copies the cell’s value to the overlying point. A terrain profile tool plots a graph with the cell’s value (elevation) on the Y axis and the distance along the section on the X axis.

Point Sampling Tool

DEM’s are often used to then update the elevation values of overlying points, for example I have used data from DEM’s to update the elevation values of address points and utilities. This isn’t as accurate as surveying each point, but it is a lot quicker! This process is also referred to image extraction, raster/vector conversion.

For this tutorial, you will need:-

  • The Point Sampling tool in QGIS is an optional plugin. You can download it by using the menus to select Plugins, Fetch Python Plugins.
  • Nasa’s srtm data, which you can download from here: http://srtm.csi.cgiar.org/
  • Some point data. If you can’t think of any, then they’re easy to create, for example use the Open Layers plugin to load Open Streetmap or Google Maps of your area, and then create points over a few cities.

I’m going to add the elevation value from the srtm rasters to a selection of UK towns and cities:-

Raster Data Extraction - UK srtm

  1. Use the menus to select Plugins, Analyses, Point Sampling Tool
  2. The Point Sampling Tool dialogue box opens. Select:-
  • The layer that contains the points to be sampled
  • The layer(s) with the field(s)/band(s) to get values from
  • The output (results) file
  • Press OK

Raster Data Extraction - Point sampling tool

The results file just contains the elevations:-

Raster Data Extraction - Elevations

It is possible to add these to the original layer:-

  • Create a buffer around the new points
  • Use the menus to select Vector, Data Management Tools, Join Attributes By location
  • Select the original points as the target and the buffer as the join layer

Another option is to update the x and y co-ordinates for both points using the Field calculator and then to match the rows in Excel on the co-ordinate column.


FOSSGIS 2013: Performance optimised wms services with QGIS server

Performance is usually a top priority for a WMS service. A recent talk at the FOSSGIS (held by Sourcepole) shows what a WMS administrator can do to optimise QGIS server performance. Finally, the performance of QGIS server is compared with UMN mapserver in two production scenarios.

Slides from FOSSGIS 2013 in Rapperswil (in german).

FOSSGIS 2013: Performance optimised wms services with QGIS server

Performance is usually a top priority for a WMS service. A recent talk at the FOSSGIS (held by Sourcepole) shows what a WMS administrator can do to optimise QGIS server performance. Finally, the performance of QGIS server is compared with UMN mapserver in two production scenarios.

Slides from FOSSGIS 2013 in Rapperswil (in german).

FOSSGIS 2013: Qualitätssicherung von Geodaten auf der Basis von Web Processing Services (WPS)

Bevor neue Daten in eine Geodaten-Infrastruktur überführt werden, müssen sie Qualitätssicherungs-Verfahren durchlaufen. Diese Prozesse können sowohl generisch, als auch speziell für ein bestimmtes Datum konfektioniert sein. Um Sicher zu stellen, dass diese Prozesse langlebig und verfügbar sind, bietet es sich an OGC-Konforme Services zu diesem Zweck zu nutzen. Deshalb werden im Vortrag die Vor- und Nachteile OGC-Konformer Web Processing Services (WPS) als Basis der Qualitätssicherung von Geodaten diskutiert.

Präsentation

FOSSGIS 2013: Qualitätssicherung von Geodaten auf der Basis von Web Processing Services (WPS)

Bevor neue Daten in eine Geodaten-Infrastruktur überführt werden, müssen sie Qualitätssicherungs-Verfahren durchlaufen. Diese Prozesse können sowohl generisch, als auch speziell für ein bestimmtes Datum konfektioniert sein. Um Sicher zu stellen, dass diese Prozesse langlebig und verfügbar sind, bietet es sich an OGC-Konforme Services zu diesem Zweck zu nutzen. Deshalb werden im Vortrag die Vor- und Nachteile OGC-Konformer Web Processing Services (WPS) als Basis der Qualitätssicherung von Geodaten diskutiert.

Präsentation

TimeManager for QGIS 2.0

As I’m sure you have already heard, QGIS 2.0 will come with a new Python API including many improvements and a generally more pythonic way of doing things. But of course that comes with a price: Old plugins (pre 2.0) won’t work anymore unless they are updated to the new version. Therefore all plugin developers are encouraged to take the time and update their plugins. An overview of changes and howto for updates is available on the QGIS wiki.

TimeManager for QGIS 2.0 will be available from day 1 of the new release. I’ve tested the usual work flows but don’t hesitate to let me know if you find any problems. The whole update process took two to three hours … sooo many signals to update … but all in all, it was far less painful than expected, thanks to everyone who contributed to the wiki update instructions!


FOSSGIS 2013: Mapfish Appserver

Mapfish Appserver is a platform for building web mapping applications using OGC standards and the Mapfish REST protocol.

Slides from FOSSGIS 2013 in Rapperswil (in german).

  • <<
  • Page 85 of 142 ( 2821 posts )
  • >>

Back to Top

Sustaining Members