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.