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 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 (DOS), je nutné jej definovat parametrem encoding.

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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/usr/bin/env python3

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 is None or prvek.attrs['psc'] != psc:
        continue
    obce_psc.add(prvek.id)
    print("{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.read_area_ids():
            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("{0}: {1}".format(prvek.attrs['psc'], prvek.attrs['nazev']))

obce.close()

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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#!/usr/bin/env python3

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 is None:
        continue
    if prvek.attrs['psc'] == psc:
        if obec_id is None:
            obec_id = prvek.id
            
        for b in prvek.boundaries():
            for n in b.read_area_ids():
                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.centroid(), attrs=(prvek.attrs['nazev'], prvek.attrs['psc']))

vystup.table.conn.commit()

vystup.close()
obce.close()

Skript ke stažení zde.

../_images/obce_psc.png

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.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
#!/usr/bin/env python3

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.