Tento víkend jsem se účastnil jako mentor akce DanubeHack a měl jsem tam možnost řešit celou řadu zajímavých úloh. Někdy se snad dostanu ke svojí malé aplikaci, která vyroběla api postavené nad Djangem a RedHatím OpenShiftem. Dneska bych se s vámi rád podělil o to, jak jsme převáděli rastrovou mapu staženou z WMS na vektory pomocí GRASSu.
Půdní mapu si vyžádal jeden z týmů a našli dokonce i funkční WMS z Pôdneho portálu Slovenska.
Půdní mapa Slovenska. Zdroj: http://www.podnemapy.sk/portal/prave_menu/podna_mapa/podna_mapa.aspx
Nejdřív jsem si založil GRASS LOCATION a MAPSET a nastavil na hranice Slovenska. Potom jsem stáhnul data z WMS služby – použil jsem formát GeoTIFF, protože JPEG díky ztrátové kompresy potrhá barvy v pixelech a místo rozumných cca 20 kategorií by jich tam bylo asi 3000. Samozřejmě jsem nastavil použitelné rozlišení regionu na 50m (původní mapa je 1:400~000):
GRASS> g.region -p projection: 99 (Krovak) zone: 0 datum: hermannskogel ellipsoid: bessel north: -1132700 south: -1334770 west: -591450 east: -165430 nsres: 50.00494927 ewres: 50.00234742 rows: 4041 cols: 8520 cells: 34429320
a stáhnu data pomocí r.in.wms
GRASS> r.in.wms url="http://sscri.vupop.sk/arcgis/services/vupop_wms/MapServer/WMSServer" wms_version="1.3.0" srs=102067 format=tiff layers=0 out=podnamapa --o --v
Jak můžete vidět, mapa je ošklivá – plná kartografických šraf a nápisů – nic, co byste mohli převézt na rozumně vypadající vektory. Přemýšlel jsem nad vhodným filtrem a nakonec z toho vypadl modul r.neighbors, který prožene přes raster okno veliké MxM pixelů, jako operátor jsem použil “nejčastěji se vyskytující hodnotu” a velikost okna zkoušel, až jsem došel k 9×9:
GRASS> r.neighbors in=podnamapa out=podnamapa_clean method=mode size=9 --o
Tím jsem se zbavil většiny šraf a nápisů (hranice se při tomhle měřítku taky trochu šouply, ale při tomhle měřítku a času, kdo by to řešil?) Pár jich ale zbylo. Při vypsání počtu pixelů za jednotlivé kategorie jsem prostě odhadnul, že ta malá čísla jsou ty špatné hodnoty – pixely, které tam nemají co dělat a které neprošly filtrem:
GRASS> r.report map=podnamapa_clean units=k +-----------------------------------------------------------------------------+ | RASTER MAP CATEGORY REPORT | |LOCATION: sk Mon Oct 19 22:30:41 2015| |-----------------------------------------------------------------------------| | north: -1132700 east: -165430 | |REGION south: -1334770 west: -591450 | | res: 50.00494927 res: 50.00234742 | |-----------------------------------------------------------------------------| |MASK: hranice in danubehack, categories 1 | |-----------------------------------------------------------------------------| |MAP: 9x9 neighborhood: mode of podnamapa2 (podnamapa2_clean in danubehack) | |-----------------------------------------------------------------------------| | Category Information | square | | #|description |kilometers| |-----------------------------------------------------------------------------| | 0| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 5.753| | 6342| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 0.010| | 8126| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 739.818| | 9664| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 0.615| |10847| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 3828.899| |13113| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 49.410| |13741| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 0.690| |14888| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 0.050| |16405| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 0.105| |17056| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 0.020| |18105| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .|10,297.153| |20250| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .|11,647.305| |20312| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 2813.458| |24319| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 901.181| |24570| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 4829.897| |25427| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 2.243| |26005| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 68.512| |26047| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 1369.785| |27482| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 3191.988| |32512| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 168.625| |32762| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 5694.238| |32767| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .| 3409.603| | *|no data. . . . . . . . . . . . . . . . . . . . . . . . . . .|37,066.504| |-----------------------------------------------------------------------------| |TOTAL |86,085.861| +-----------------------------------------------------------------------------+
Pomocí r.reclass jsem se jich zbavil:
GRASS> r.reclass in=podnamapa_clean out=podnamapa_clean2 rules=- << EOF 0 = NULL 6342 = NULL 8126 = 8126 9664 = NULL 10847 = 10847 13113 = 13113 13741 = NULL 14888 = NULL 16405 = NULL 17056 = NULL 18105 = 18105 20250 = 20250 20312 = 20312 24319 = 24319 24570 = 24570 25427 = 25427 26005 = 26005 26047 = 26047 27482 = 27482 32512 = 32512 32762 = 32762 32767 = 3767 * = NULL EOF
Převod raster na vektor je už hračka
GRASS> r.to.vect in=podnamapa_clean2 out=podnamapa type=area
a aby nebyly na polygonech oškolivé pixely, tak jsem to ještě prohnal generalizací
GRASS> v.generalize in=podnamap out=podnamapa_generalize --o method=douglas threshold=100
Výsledná mapa nemá kategorie, to si ale už účastníci Hackathonu doplnili z původního datového zdroje ručně.
A že jste to vy, vrazil jsem to ve formátu TopoJSON na GitHub, abyste ti to mohli prohlídnout a stáhnout (doplnění kategorií za domácí úkol).
Výsledek mi zabral asi tak hodinu zkoušení. Má to chyby: např. díry po té reklasifikaci, že, ale šlo o čas a “dostatečně dobré” řešení.
