Implementace skriptu pro GRASS

Skript bude v první fázi vypisovat pro dané PSČ seznam obcí a jejich sousedy dle PSČ.

Příprava dat

Vektorová mapa obce_polygon z mapsetu ruian neobsahuje v atributové tabulce informace o PSČ. Tuto informaci budeme muset nejprve doplnit. PSČ obcí ke stažení zde jako DBF tabulka.

V aktuálním mapsetu si vytvoříme kopii původní mapy (g.copy) a k její atributové tabulce připojíme tabulku s PSČ.

g.copy vector=obce_polygon@ruian,obce

Vstupní tabulku s PSČ nejprve pomocí db.in.ogr naimportujeme:

db.in.ogr input=zv_pcobc.dbf encoding=cp852

Poznámka

Modul předpokládá kódování znaků UTF-8. Vzhledem k tomu, že jsou vstupní data v kódování CP852, je nutné jej definovat parametrem encoding.

Důležité

Ve verzích GRASS 7.0.2 a nižší modul db.in.ogr parametr encoding nemá. V tomto případě definujeme kódování pomocí proměnné prostředí SHAPE_ENCODING.

SHAPE_ENCODING=cp852 db.in.ogr input=zv_pcobc.dbf

Tabulka bohužel neobsahuje kód obce, podle kterého bysme mohli data pohodlně pomocí modulu v.db.join připojit k atributové tabulce obcí, viz. výpis modulu db.select.

db.select sql="select * from zv_pcobc_dbf limit 1"

NAZCOBCE|PSC|NAZPOST|KODOKRESU|NAZOKRESU|NAZOBCE|KOD
Abertamy|36235|Abertamy|3403|Karlovy Vary|Abertamy|

Vytvoříme pomocnou tabulku pomocí db.execute, která bude obsahovat dva sloupce: kód obce a PSČ.

db.execute sql="create table obce_psc as select o.kod,z.psc from obce as o \
 join zv_pcobc_dbf as z on o.okreskod = z.KODOKRESU and o.nazev = z.NAZOBCE"

Poznámka

Tento postup je nutný, neboť databáze SQLite (která je pro systém GRASS vychozí při ukládání atributových dat, viz školení pro začátečníky), nepodporuje JOIN v příkazu UPDATE. Pokud používate namísto toho např. PostgreSQL, tak se operace zjednoduší.

Takto vytvořenou tabulku již lze k atributové tabulce vektorové mapy obcí pomocí v.db.join připojit.

v.db.join map=obce column=kod other_table=obce_psc other_column=kod

Výsledek si můžeme zkontrolovat pomocí modulu v.db.select.

v.db.select map=obce column=kod,nazev,psc

...
532401|Jarpice|27371
532428|Jemníky|27401
548103|Kámen|39413
...

Implementace testovací verze skriptu

Podrobnější informace v kapitole Přístup k vektorovým datům.

#!/usr/bin/env python

from grass.pygrass.vector import VectorTopo

psc = '41115'

obce = VectorTopo('obce')
obce.open('r')

print ("Seznam obci s PSC {}:".format(psc))
obce_psc = set()
for prvek in obce.viter('areas'):
    if prvek.attrs['psc'] != psc:
        continue
    obce_psc.add(prvek.id)
    print (u"{0}: {1}".format(psc, prvek.attrs['nazev']))
    
sousede = set()
for prvek in obce.viter('areas'):
    if prvek.id not in obce_psc:
        continue

    for b in prvek.boundaries():
        for n in b.get_left_right():
            if n != -1 and n != prvek.id:
               sousede.add(n)

print ("Seznam sousednich obce:")
for prvek in obce.viter('areas'):
    if prvek.id not in sousede or \
       prvek.attrs['psc'] == psc:
        continue
    
    print (u"{0}: {1}".format(prvek.attrs['psc'], prvek.attrs['nazev']))

obce.close()

Varování

Mezi verzemi GRASS 7.1 a GRASS 7.0 se API PyGRASS částečně změnilo. Místo funkce get_left_right() (24) použijte v GRASS 7.1 read_area_ids().

Poznámka

Skript je koncipován jako ukázka, uvedená implementace má vážné nedostatky, např. nepoužívá vůbec prostorový index. Místo toho všechny prvky prochází sekvenčně.

Skript ke stažení zde.

Výstup může vypadat následovně:

Seznam obci s PSC 41115:
41115: Děčany
41115: Podsedice
41115: Jenčice
41115: Lkáň
41115: Chodovlice
41115: Dlažkovice
41115: Sedlec
41115: Třebívlice
Seznam sousednich obce:
44001: Želkovice
41002: Velemín
41113: Třebenice
41804: Lukov
41002: Úpohlavy
41002: Slatina
41114: Vlastislav
43921: Koštice
43922: Chožov
43926: Libčeves
41002: Vchynice
41002: Čížkovice
41116: Klapý
41113: Třebenice
41757: Hrobčice

Druhá verze skriptu

Druhá verze skriptu se bude lišit od první tím, že namísto textového výpisu, vytvoří novou vektorovou mapu se spojeným polygonem obcí s definovaným PSČ a zároveň bude obsahovat sousední obce.

#!/usr/bin/env python

from grass.pygrass.vector import VectorTopo

psc = '41115'

obce = VectorTopo('obce')
obce.open('r')

vystup = VectorTopo('obce_psc_{}'.format(psc))
vystup.open('w', tab_cols=[('cat',       'INTEGER PRIMARY KEY'),
                           ('nazev',     'TEXT'),
                           ('psc',       'INTEGER')])

obec_id = None
obce_psc = set()
for prvek in obce.viter('areas'):
    if prvek.attrs['psc'] == psc:
        if obec_id is None:
            obec_id = prvek.id
            
        for b in prvek.boundaries():
            for n in b.get_left_right():
                if n != -1 and n != obec_id:
                    obce_psc.add(n)
obce_psc.add(obec_id)

hranice = list()
for prvek in obce.viter('areas'):
    if prvek.id not in obce_psc:
        continue

    for b in prvek.boundaries():
        if b.id not in hranice:
            hranice.append(b.id)
            vystup.write(b, attrs=(None, None))

    vystup.write(prvek.get_centroid(), attrs=(prvek.attrs['nazev'], prvek.attrs['psc']))

vystup.table.conn.commit()

vystup.close()
obce.close()

Varování

Mezi verzemi GRASS 7.1 a GRASS 7.0 se API PyGRASS částečně změnilo. Místo funkce get_centroid() (38) použijte v GRASS 7.1 centroid().

Skript ke stažení zde.

Obr. 15 Vizualizace výsledku pro PSČ 41115.

Poznámka

Úloha by se dala samozřejmě řešit poměrně jednoduše namísto PyGRASS kombinací standardních modulů systému GRASS v.extract a v.select.

#!/usr/bin/env python

from grass.pygrass.modules import Module

psc = '41115'

Module('v.extract', input='obce', output='obce1',
       where="psc = '{}'".format(psc))
Module('v.select', ainput='obce', binput='obce1',
       output='obce_psc_{}'.format(psc),
       operator='overlap', overwrite=True)
Module('g.remove', type='vector', name='obce1', flags='f')

Skript ke stažení zde.