Geocoding with GDAL: building an API

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
It's only fair to share...Tweet about this on TwitterShare on FacebookShare on LinkedInPin on PinterestShare on RedditShare on Google+

Trackbacks/Pingbacks

  1.  ItOpen – Open Web Solutions, WebGis Development » Blog Archive » geocoding with GDAL