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#!/usr/bin/env python3
 2
 3from grass.pygrass.vector import VectorTopo
 4
 5psc = '41115'
 6
 7obce = VectorTopo('obce')
 8obce.open('r')
 9
10print("Seznam obci s PSC {}:".format(psc))
11obce_psc = set()
12for prvek in obce.viter('areas'):
13    if prvek.attrs is None or prvek.attrs['psc'] != psc:
14        continue
15    obce_psc.add(prvek.id)
16    print("{0}: {1}".format(psc, prvek.attrs['nazev']))
17    
18sousede = set()
19for prvek in obce.viter('areas'):
20    if prvek.id not in obce_psc:
21        continue
22
23    for b in prvek.boundaries():
24        for n in b.read_area_ids():
25            if n != -1 and n != prvek.id:
26               sousede.add(n)
27
28print("Seznam sousednich obce:")
29for prvek in obce.viter('areas'):
30    if prvek.id not in sousede or \
31       prvek.attrs['psc'] == psc:
32        continue
33    
34    print("{0}: {1}".format(prvek.attrs['psc'], prvek.attrs['nazev']))
35
36obce.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#!/usr/bin/env python3
 2
 3from grass.pygrass.vector import VectorTopo
 4
 5psc = '41115'
 6
 7obce = VectorTopo('obce')
 8obce.open('r')
 9
10vystup = VectorTopo('obce_psc_{}'.format(psc))
11vystup.open('w', tab_cols=[('cat',       'INTEGER PRIMARY KEY'),
12                           ('nazev',     'TEXT'),
13                           ('psc',       'INTEGER')])
14
15obec_id = None
16obce_psc = set()
17for prvek in obce.viter('areas'):
18    if prvek.attrs is None:
19        continue
20    if prvek.attrs['psc'] == psc:
21        if obec_id is None:
22            obec_id = prvek.id
23            
24        for b in prvek.boundaries():
25            for n in b.read_area_ids():
26                if n != -1 and n != obec_id:
27                    obce_psc.add(n)
28obce_psc.add(obec_id)
29
30hranice = list()
31for prvek in obce.viter('areas'):
32    if prvek.id not in obce_psc:
33        continue
34
35    for b in prvek.boundaries():
36        if b.id not in hranice:
37            hranice.append(b.id)
38            vystup.write(b, attrs=(None, None))
39
40    vystup.write(prvek.centroid(), attrs=(prvek.attrs['nazev'], prvek.attrs['psc']))
41
42vystup.table.conn.commit()
43
44vystup.close()
45obce.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#!/usr/bin/env python3
 2
 3from grass.pygrass.modules import Module
 4
 5psc = '41115'
 6
 7Module('v.extract', input='obce', output='obce1',
 8       where="psc = '{}'".format(psc))
 9Module('v.select', ainput='obce', binput='obce1',
10       output='obce_psc_{}'.format(psc),
11       operator='overlap', overwrite=True)
12Module('g.remove', type='vector', name='obce1', flags='f')

Skript ke stažení zde.