This is the second part (read the first part) of my GDAL geocoding experiments, here I’m trying to build an API that makes it simpler to use GDAL geocoding services from python. The good ideas are proudly stolen from the GDAL binding available in the Django project, bad ideas are all mine 🙂
Setting up the bindings
In this startup code, we are importing ctypes and six (for python 2/3 compatibility) and telling ctypes that some functions are returning (char pointers) char * , this is important to get back python strings from the function calls.
Since ctypes does not support enums, a list of important GDAL constants is also provided.
#!/bin/env python
"""
Testing geocoding GDAL implementation: building an API.
Most of the ideas here are borrowed from Django GDAL
bindings.
"""
import six
from ctypes import *
lgdal = CDLL('libgdal.so')
# Setup bindings
lgdal.OGR_F_GetFieldAsString.restype = c_char_p
lgdal.OGR_Fld_GetNameRef.restype = c_char_p
lgdal.OGR_G_ExportToWkt.argtypes = [c_void_p, POINTER(c_char_p)]
# ctypes does not support enums...
# enum OGRFieldType
OFTInteger = 0
OFTIntegerList = 1
OFTReal = 2
OFTRealList = 3
OFTString = 4
OFTStringList = 5
OFTWideString = 6
OFTWideStringList = 7
OFTBinary = 8
OFTDate = 9
OFTTime = 10
OFTDateTime = 11
# enum OGRwkbGeometryType
wkbUnknown = 0
wkbPoint = 1
wkbLineString = 2
wkbPolygon = 3
wkbMultiPoint = 4
wkbMultiLineString = 5
wkbMultiPolygon = 6
wkbGeometryCollection = 7
wkbNone = 100
wkbLinearRing = 101
wkbPoint25D = 0x80000001
wkbLineString25D = 0x80000002
wkbPolygon25D = 0x80000003
wkbMultiPoint25D = 0x80000004
wkbMultiLineString25D = 0x80000005
wkbMultiPolygon25D = 0x80000006
wkbGeometryCollection25D = 0x80000007
Building the API
In this part of the script we first create a base class to hold the c pointer to the underlying GDAL object and then we create the wrapper classes. The bindings are not complete because not only the process would be longer and boring but also the resulting API would have been harder to use and less “pythonic”.
class GDALBase(object):
"""
Base object for GDAL objects that has a pointer access property
that controls access to the underlying C pointer.
"""
# Initially the pointer is NULL.
_ptr = None
# Default allowed pointer type.
ptr_type = c_void_p
# Pointer access property.
def _get_ptr(self):
# Raise an exception if the pointer isn't valid don't
# want to be passing NULL pointers to routines --
# that's very bad.
if self._ptr: return self._ptr
else: raise GDALException('GDAL %s pointer no longer valid.' \
% self.__class__.__name__)
def _set_ptr(self, ptr):
# Only allow the pointer to be set with pointers of the
# compatible type or None (NULL).
if isinstance(ptr, six.integer_types):
self._ptr = self.ptr_type(ptr)
elif ptr is None or isinstance(ptr, self.ptr_type):
self._ptr = ptr
else:
raise TypeError('Incompatible pointer type')
ptr = property(_get_ptr, _set_ptr)
class OGRGeocodingSession(GDALBase):
"""The geocoding session: this holds the configuration"""
def __init__(self, *args):
session_config_type = c_char_p * len(args) # Holds config strings
session_config = session_config_type(*args)
self.ptr = lgdal.OGRGeocodeCreateSession(byref(session_config))
def __del__(self):
"""Destroy session"""
lgdal.OGRGeocodeDestroySession(self.ptr)
Here comes the main wrapper class, this class holds a reference to the session so that we can hide this complexity from the user and forget about destroying the session and layer pointers when they are not needed (see __del__ functions).
class OGRGeocoder(GDALBase):
"""The geocoder wrapper: holds session and layer instances"""
_session = None
_layer = None
def __init__(self, *args):
"""Configuration options are available in args"""
self._session = OGRGeocodingSession(*args)
if not self._session:
raise GDALException('GDAL cannot create geocoder session.')
def __del__(self):
"""Destroy session and free layer"""
if self._layer is not None:
lgdal.OGRGeocodeFreeResult(self._layer)
# Session is deleted automatically whenever it goes out of scope
def _layer_to_python(self):
"""Converts GDAL geocoding results into a python iterator which
returns its features as a dictionary"""
if not self._layer:
return []
else:
def _result_row(layer, num_results, fields):
"""Fetch row"""
num = 0
while num < num_results:
feature = lgdal.OGR_L_GetNextFeature(layer)
num += 1
row = {'__num': num}
# Scan all fields
idx = 0
for f in fields:
fname = lgdal.OGR_Fld_GetNameRef(f)
ftype = lgdal.OGR_Fld_GetType(f)
if ftype == OFTInteger:
value = lgdal.OGR_F_GetFieldAsInteger(feature, idx)
elif ftype == OFTReal:
value = lgdal.OGR_F_GetFieldAsDouble(feature, idx)
else:
value = lgdal.OGR_F_GetFieldAsString(feature, idx)
row[fname] = value
idx += 1
# Adds geometry
geometry = lgdal.OGR_F_GetGeometryRef(feature)
if geometry is not None:
buf = c_char_p()
lgdal.OGR_G_ExportToWkt(geometry, buf)
row['geometry'] = buf.value
yield row
num_results = lgdal.OGR_L_GetFeatureCount(self._layer)
field_definitions = lgdal.OGR_L_GetLayerDefn(self._layer)
num_fields = lgdal.OGR_FD_GetFieldCount(field_definitions)
fields = [lgdal.OGR_FD_GetFieldDefn(field_definitions, i) \
for i in range(num_fields)]
return _result_row(self._layer, num_results, fields)
def geocode(self, address):
"""Calls the geocoding service with the given search string"""
if self._layer is not None:
lgdal.OGRGeocodeFreeResult(self._layer)
self._layer = lgdal.OGRGeocode(self._session.ptr, address, None, None)
return self._layer_to_python()
def reverse(self, lon, lat):
"""Calls the reverse geocoding service with the given lon/lat"""
if self._layer is not None:
lgdal.OGRGeocodeFreeResult(self._layer)
self._layer = lgdal.OGRGeocodeReverse(self._session.ptr, c_double(lon), \
c_double(lat), None)
return self._layer_to_python()
Finally, we are going to test our classes:
# Test it!
if __name__ == "__main__":
def _print_results(results, title):
# Formatted print
try:
for result in results:
print(('=' * 10) + title + ('=' * 20))
for k in result:
print("%s\t%s" % (k, result[k]))
except StopIteration:
print("No results for %s!" % title)
gc = OGRGeocoder()
results = gc.geocode('Milano')
_print_results(results, 'Milano')
results = gc.geocode('Merano')
_print_results(results, 'Merano')
# Test a not existing location
results = gc.geocode('')
_print_results(results, '')
# Test reverse
results = gc.reverse(11.149171, 46.6731974)
_print_results(results, 'Reverse')
The results (extract) follows:
==========Milano====================
place_id 97562386
hamlet
county MI
postcode
country_code it
village
city Milano
display_name Milano, MI, LOM, Italia
locality
lon 2
state LOM
boundingbox 45.3867378234863,45.5358505249023,9.04065799713135,9.27803993225098
type city
importance 0.87585157580298
__num 1
lat 2
class place
icon http://nominatim.openstreetmap.org/images/mapicons/poi_place_city.p.20.png
geometry POLYGON ((9.040658 45.4474642,9.0411095 45.4477548,9.0439066 45.4486802, .....))
country Italia
osm_id 44915
osm_type relation
place_rank 16
==========Milano====================
place_id 98024954
hamlet
county Provincia de Salamanca
postcode
country_code es
village
city El Milano
display_name El Milano, Provincia de Salamanca, Castilla y León, España
locality
lon 2
state Castilla y León
boundingbox 41.0766944885254,41.1206321716309,-6.64893770217896,-6.56636381149292
type administrative
importance 0.53538498085989
__num 2
lat 2
class boundary
icon http://nominatim.openstreetmap.org/images/mapicons/poi_boundary_administrative.p.20.png
geometry POLYGON ((-6.6489375 41.1054302,-6.6450082 41.1062842,-6.6377713 41.1086927, .....))
country España
osm_id 344499
osm_type relation
place_rank 16
==========Merano====================
house
place_id 97168116
county BZ
city_district
postcode
country_code it
osm_id 47283
city Meran - Merano
display_name Meran - Merano, BZ, TAA, Italia
lon 2
state TAA
station
boundingbox 46.6168327331543,46.6944541931152,11.1364288330078,11.2397155761719
type administrative
state_district
importance 0.65497157559577
hotel
__num 1
suburb
lat 2
class boundary
icon http://nominatim.openstreetmap.org/images/mapicons/poi_boundary_administrative.p.20.png
bar
geometry POLYGON ((11.1364291 46.6896835,11.1372572 46.6899335,11.1375407 46.6900137,....))
country Italia
osm_type relation
place_rank 16
address29
road
==========Merano====================
house
place_id 16545865
county BZ
city_district
postcode 39012
country_code it
osm_id 1455273173
city Meran - Merano
display_name Meran/Merano, Meran Bahnhof, Meran - Merano, BZ, TAA, 39012, Italia
lon 2
state TAA
station Meran/Merano
boundingbox 46.6731948852539,46.6731986999512,11.1491708755493,11.1491718292236
type station
state_district
importance 0.24748918685134
hotel
__num 2
suburb
lat 2
class railway
icon http://nominatim.openstreetmap.org/images/mapicons/transport_train_station2.p.20.png
bar
geometry POINT (11.149171 46.6731974)
country Italia
osm_type node
place_rank 30
address29
road Meran Bahnhof
==========Reverse====================
town Meran - Merano
display_name 8, Bahnhofsplatz / Piazza Stazione, Meran - Merano, BZ, TAA, 39012, Italia
house_number 8
country Italia
place_id 9172391483
lon 2
state TAA
__num 1
county BZ
geometry POINT (11.1492877 46.6732366)
osm_type node
postcode 39012
country_code it
osm_id 2755814453
lat 2
road Bahnhofsplatz / Piazza Stazione
This is the second part (read the first part) of my GDAL geocoding experiments, here I’m trying to build an API that makes it simpler to use GDAL geocoding services from python. The good ideas are proudly stolen from the GDAL binding available in the Django project, bad ideas are all mine 🙂
