Geokódování

Tato stránka shrnuje poznámky k procesu geokódování, pomocí kterého se textovým adresním záznamům v tabulce/databázi přiřadí prostorové souřadnice X,Y s využitím uliční sítě.

Rozlišujeme:

• dopředné geokódování (forward geocoding), kdy na základě adres míst určujeme prostorové souřadnice
• zpětné geokódování (backward geocoding), kdy naopak na základě lokalizace geoprvků odvozujeme např. jejich adresy v uliční síti

Existuje několik geokódovacích služeb. Mezi nejrozšířenější patří Nominatim, který je založen na datech z projektu OpenStreetMap.

Příklad skriptu pro Nominatim

Níže uvedený skript na základě adres uložených v databázi odvodí souřadnice jednotlivých prvků v souřadnicovém systému WGS-84 (EPSG 4326) a uloží je do dané tabulky (sloupce lat, long). Kvalitu geokódování uloží do sloupce level jako hodnoty:

• 0 - geokódováno na ulici, čp, část obce a obec
• 1 - geokódováno na ulici část obce a obec
• 2 - geokódováno na část obce a obec
• 3 - geokódováno na obec
• 4 - geokódováno na obec (první část před čárkou)

Příklad vstupních dat v DB:

 id1 |    ulice    | cp_co |  cast_obce  |      obec      |  psc  
-----+-------------+-------+-------------+----------------+-------
   1 | Rákosníkova | 162   | Neveklov    | Neveklov       | 25756
   2 | Lesní cesta | 236/5 | Klíše       | Ústí nad Labem | 40001
   3 | Svojsíkova  | 833   | Český Těšín | Český Těšín 1  | 73701
   4 | Vídeňská    | 9     | Štýřice     | Brno 39        | 63900
   5 | Hradská-KD  | 723/8 | Loštice     | Loštice        | 78983
...

Spuštění skriptu:

python chess.py
1 -> 49.7541028,14.5363398 (0)
2 -> 50.6730414,14.0141232 (0)
3 -> 49.74285935,18.6228764915793 (3)
4 -> 49.1817338,16.5951966 (1)
5 -> 49.7449432,16.9281353 (3)
...

Skript vypisuje na standardní výstup informaci o id -> lat,long (level). Výsledek je SQL dávka, kterou můžeme nahrát do databáze.

V našem případě může výsledná dávka vypadat následovně:

BEGIN;
update tabulka set lat=49.7541028,lon=14.5363398,level=0 where id1=1;
update tabulka set lat=50.6730414,lon=14.0141232,level=0 where id1=2;
update tabulka set lat=49.7428593,lon=18.6228764,level=3 where id1=3;
update tabulka set lat=49.1817338,lon=16.5951966,level=1 where id1=4;
update tabulka set lat=49.7449432,lon=16.9281353,level=3 where id1=5;
...
COMMIT;

Výsledek nahrajeme do databáze:

psql -d databaze -f davka.sql

Příklad výsledku geokódování:

 id1 |    ulice    | cp_co |  cast_obce  |      obec      |  psc  |     lat     |       lon        | level 
-----+-------------+-------+-------------+----------------+-------+-------------+------------------+-------
   1 | Rákosníkova | 162   | Neveklov    | Neveklov       | 25756 |  49.7541028 |       14.5363398 |     0
   2 | Lesní cesta | 236/5 | Klíše       | Ústí nad Labem | 40001 |  50.6730414 |       14.0141232 |     0
   3 | Svojsíkova  | 833   | Český Těšín | Český Těšín 1  | 73701 |  49.7428593 |       18.6228764 |     3
   4 | Vídeňská    | 9     | Štýřice     | Brno 39        | 63900 |  49.1817338 |       16.5951966 |     1
   5 | Hradská-KD  | 723/8 | Loštice     | Loštice        | 78983 |  49.7449432 |       16.9281353 |     3

Ukázka skriptu (Python)

import psycopg2
import urllib2
import json
import time
import sys
import os

# nazev databaze
DB='databaze'
# sloupce v tabulce obsahujici klic, ulici, c.p., cast obce a obec
COLUMNS = ['id1', 'ulice', 'cp_co', 'cast_obce', 'obec']
# nazev tabulky
TABLE='tabulka'
# cesta k vystupnimu souboru
OUTPUT='davka.sql'

def get_data(db, table, columns):
    conn = psycopg2.connect("dbname=%s" % db)
    cur = conn.cursor()

    cur.execute("SELECT %s FROM %s WHERE %s IS NOT NULL ORDER BY %s" % (','.join(columns), table,
                                                                        columns[-1], columns[0]))
    data = []
    for rec in cur.fetchall():
        row = {}
        for idx in range(len(columns)):
            row[columns[idx]] = rec[idx]
        data.append(row)
        
    cur.close()
    conn.close()
    
    return data

def geocode(data, key, output, ocols=('lat', 'lon', 'level')):
    fd = open(output, 'w')
    fd.write('BEGIN;%s' % os.linesep)
    for addr in data:
        level=0
        while True:
            try:
                data = json.loads(try_it(addr, level))[0]
            except:
                level += 1
            break

        lat = data['lat']
        lon = data['lon']
        
        sys.stdout.write("%s -> %s,%s (%s)%s" % (addr[key], lat, lon, level, os.linesep))
        fd.write('update chess_oddil2 set %s=%s,%s=%s,%s=%d where %s=%s;%s' % \
                 (ocols[0], lat, ocols[1], lon, ocols[2], level, key, addr[key], os.linesep))
        
    fd.write('COMMIT;%s' % os.linesep)
    fd.close()
        
def try_it(addr, level=0):
    qstring = ''
    if level < 2 and addr['ulice']:
        qstring += addr['ulice'].replace(' ', '+')
        if level == 0 and addr['cp_co']:
            if '/' in addr['cp_co']:
                qstring += '+' + addr['cp_co'].split('/')[0]
            else:
                qstring += '+' + addr['cp_co']
        qstring += '+' + addr['obec'].replace(' ', '+')
    else:
        if level < 3:
            qstring += addr['obec'].replace(' ', '+')
        else:
            if ',' in addr['obec']:
                qstring += addr['obec'].split(',')[0].replace(' ', '+')
            else:
                qstring += addr['obec'].split(' ', 1)[0]
    url = 'http://nominatim.openstreetmap.org/search?q=' + qstring + '&format=json&limit=1'
    
    response = urllib2.urlopen(url)
    jsondata = response.read()
    
    return jsondata
    
if __name__ == "__main__":
    geocode(data=get_data(DB, TABLE, COLUMNS), key=COLUMNS[0], output=OUTPUT)
    sys.exit(0)