This post shows how we can query geographic data using a simple Python client pointing to a GIS web service: SnowFloat.

The data file here is the Vermont towns geographic data file. Its name is “Boundary_BNDHASH_region_towns”. We first upload the data file to the web service. The web service here is a layer on top of GeoDjango on top of PostGIS.

Layer

Let’s retrieve the layer object:

>>> layer = client.get_layers(name='Boundary_BNDHASH_region_towns')[0]

Number of features:

>>> layer.num_features
255

Layer’s fields:

>>> layer.fields
[{u'name': u'gdb_vcgi_v', u'type': u'real'},
 {u'name': u'perimeter', u'type': u'real'},
 {u'name': u'towns_', u'type': u'integer'},
 {u'name': u'towns_id', u'type': u'integer'},
 {u'name': u'fips6', u'type': u'integer'},
 {u'name': u'townname', u'size': 18, u'type': u'string'},
 {u'name': u'townnamemc', u'size': 18, u'type': u'string'},
 {u'name': u'cnty', u'type': u'integer'},
 {u'name': u'shape_area', u'type': u'real'},
 {u'name': u'shape_len', u'type': u'real'}]

Layer’s extent:

>>> layer.extent
[424788.84, 25211.8370313, 581562.0, 279799.0]

Features

Towns starting with the letter ‘B':

>>> features = layer.get_features(field_townnamemc_startswith='B')
>>> len(features)
30
>>> features[0].fields['townnamemc']
Berkshire

Towns less than 10 kms from the town of Burlington:

>>> burlington = layer.get_features(field_townnamemc_exact='Burlington')
>>> features = layer.get_features(query='distance_lte',
...                               geometry=burlington,
...                               distance=10000)
>>> len(features)
12
>>> [feature.fields['townnamemc'] for feature in features]
[u'Milton',
 u'South Hero',
 u'St. George',
 u'Charlotte',
 u'Colchester',
 u'Winooski',
 u'Shelburne',
 u'Hinesburg',
 u'Burlington',
 u'Essex',
 u'South Burlington',
 u'Williston']

Ten furthest towns + their distance to Burlington:

>>> features = layer.get_features(spatial_operation='distance',
...                               spatial_geometry=burlington.geometry,
...                               order_by=('spatial',),
...                               query_slice=(0, 10))
>>> ['%s - %d kms' % (feature.fields['townnamemc'], feature.spatial / 1000) for feature in features]
[u'Vernon - 187 kms',
 u'Guilford - 184 kms',
 u'Halifax - 182 kms',
 u'Whitingham - 180 kms',
 u'Stamford - 179 kms',
 u'Pownal - 178 kms',
 u'Brattleboro - 177 kms',
 u'Readsboro - 177 kms',
 u'Marlboro - 172 kms',
 u'Dummerston - 170 kms']

Towns inside the Robert T. Stafford White Rocks area:

>>> features = layer.get_features(query='intersects', geometry=forest_feature.geometry)
>>> [feature.fields['townnamemc'] for feature in features]
[u'Wallingford', u'Mount Holly', u'Mount Tabor', u'Weston', u'Dorset', u'Peru']

As you can see, the goal here is to build a light tool for easy access to geographic data. As I am continuing to work on the web service back-end, I would really appreciate getting feedback from people working in the GIS world.