Feature Query Howto

Version: 0.7
Author: Sean Gillies
Revision: 1.1
Date: 2005-02-15

Introduction

Version 0.7 of the Python Cartographic Library introduces a feature querying interface that is inspired by the OGC's WFS GetFeature specification. Unlike MapServer's mapscript, with PCL we query data stores directly (depicted below)

  [ Results ]
       ^
       |
---------------
 Query Filters
---------------
       ^
       |
 +-----------+
 | DataStore |
 +-----------+

rather than querying a layer or map. In this release, we have querying for local disk data stores only.

Note

the new querying features require the Python modules from the GDAL project. Please see the 0.7 release notes on the project download site or the PCL/RELEASE.txt file included in version 0.7 for further directions regarding GDAL.

Before the 1.0 release we will develop querying for external data stores on an as-needed basis. I presume that PostGIS and WFS data stores will be the first.

Query and Results

Forming a query

A feature query is mediated by the data.Query class. Instances of Query specify by name the desired feature properties, the spatial and attribute filters, and the environment for evaluating filter expressions. See the API documentation

http://zmapserver.sourceforge.net/PCL/api/index.html

for details.

Query results

The query object is passed to the features method of a data store which returns a list of feature results:

>>> q = data.Query()
>>> results = datastore.features(q)

or to the store's iterfeatures method which returns an iterator over query results:

>>> q = data.Query()
>>> for result in datastore.iterfeatures(q):
>>>     ...

Filters

Python attribute and spatial filters

Filtering of features is specified by the filter.Filter class and its subclasses. The important subclass for the 0.7 release is PythonFilter for local data stores. This filter uses Python expressions and evaluates them within a local context which includes the feature geometry and properties. The feature properties can be referenced as attributes of an object named "f". For example, the following PythonFilter expression:

'f.FIPS_CNTRY in ["US", "GM", "FR", "FI"]'

could be used with the world borders dataset (used in the ZCO manual) from

http://mappinghacks.com/data/

to select features from four countries that are home to ZMapServer and ZCO developers. To select features from countries that have a population density greater than 300 persons per square kilometer, you could use:

'f.POP_CNTRY/f.AREA > 300'

The feature geometry itself is known within the PythonFilter expression as "g". In 0.7 this object has two methods that can be used in expressions: BBOX(), and intersects(). The BBOX method can be used when querying local disk data stores to efficiently find features whose geometries do relate to the specified bounding box. By "relate" I mean within or intersecting with the box. Arguments to BBOX must be four literal floating point values in the standard minx, miny, maxx, maxy order. The expression:

'g.BBOX(-20.0,25.0,60.0,75.0)'

can be used to select features from a few Mediterranean and European countries.

Note

a bounding box query expression will always be evaluated prior to any other attribute or spatial expressions and so can be used to improve the performance of a query.

To exploit the intersects method and the rest of the upcoming spatial predicates, you will need to modify the filter expression environment or namespace:

>>> from cartography import data, spatial
>>> from cartography.filter import PythonFilter
>>> ds = data.DiskFeatureStore('/var/wms_data/world/world_borders.shp')
>>> f = PythonFilter('g.intersects(spatial.Point(10.0, 40.0))')
>>> q = data.Query(filter=f, spatial=spatial)
>>> ds.features(q)
[(<cartography.spatial.Polygon instance at 0xf703394c>, [111, 'IT', 'Italy', 301230.0, 58057477.0])]

In the example above, we add the spatial module to the filter evaluation environment by passing it in through the Query class's params mapping. Until we have done this, the only names in the filter's namespace are the standard Python built-ins and our "f" and "g". We can also create the point outside of the expression to improve performance:

>>> f = PythonFilter('g.intersects(p)')
>>> query_point = spatial.Point(10.0, 40.0)
>>> q = data.Query(filter=f, p=query_point)
>>> ds.features(q)
[(<cartography.spatial.Polygon instance at 0xf703394c>, [111, 'IT', 'Italy', 301230.0, 58057477.0])]

Non-Python filters and external data stores

The examples above apply only to local disk data stores. I intend for the PCL project to develop a more general set of classes for creating Filters than can be used with external data stores. This would place all the filtering responsibility in the hands of the external service, and our PCL classes would simply (!) serialize themselves into XML for querying an external WFS or SQL for querying a PostGIS (or other) database.

I do not intend to work on munging PythonFilter expressions into XML or SQL.