Denna vecka tänkte jag titta lite på hur man kan utöka funktionaliteten i sitt GIS program med Pythonskript och extra programbibliotek.
Det kan ju vara så att det GIS program man använder saknar en funktionalitet som man behöver av olika skäl, men att det finns exempelvis terminalprogram som skulle kunna lösa det man vill göra, fast dessa är ju inte alltid så användarvänliga.
I tre delar denna veckan så går jag igenom ett exempel på hur man kan använda ett sådant terminalprogram för att bygga in funktionaliteten i sitt GIS.
Idag blir det grunder och kod, sedan integrerar jag skriptet i QGIS och avslutar med ArcMap.
Är du utvecklare av Pythonskript för GIS och orkar hålla ut till på onsdag så kommer jag att bjuda på något som kanske är en nyhet för dig, och som kommer att göra livet mycket trevligare när det gäller exempelvis att öka kundbasen.
Allt började (den här gången) med att jag inte hittade ett verktyg för att dela upp en stor rasterbild i mindre, angränsande, lika stora rasterbilder. Däremot så fanns det många som hänvisade till ett GDAL kommando som fungerar såhär:
Ja, självklart!!!
Börjar man bryta i kommandot så ser man ganska snabbt att det är ett kommando som tar minst två argument. Ett in-dataset och ett ut-dataset. Men granskar man lite närmare så ser man också att om det jag vill skall göras så måste man ge kommandot en gång (med massor av parametrar) för varje liten rasterbild som skall bli resultatet.
Detta gör att det lämpar sig mycket väl för att skapa ett skript.
Skriptet skall ta tre argument:
- Rasterfilen som skall skäras upp
- Hur stora skall ut-datafilerna vara (pixlar i sida)
- Var skall ut-datafiler sparas (katalog)
Det fanns flera förslag på hur detta skrivs i Python och jag valde ett av dessa, som jag sedan anpassade för att det skulle passa mig bättre.
#!/usr/bin/env python import os, sys from osgeo import gdal r_src = gdal.Open(sys.argv[1]) # läs in det första argumentet som rasterfil i variabeln 'r_src' dst_size = int(sys.argv[2]) # hur stora skall resultatfilerna vara (pixlar i sida) r_dst = sys.argv[3] # katalog där de resulterande filerna sparas width = r_src.RasterXSize # placera rasterbildens pixelbredd i variabeln 'width' height = r_src.RasterYSize # placera rasterbildens pixelhöjd i variabeln 'height' print width, 'x', height # skriv ut bredd och höjd till konsolen for i in range(0, width, dst_size): # upprepa det antal gånger som 'dst_size' ryms i bildbredden (anges med 'i') for j in range(0, height, dst_size): # upprepa det antal gånger som 'dst_size' ryms i bildhöjden (anges med 'j') w = min(i+dst_size, width) - i # tar det som är minst i+dst_size eller rasterbredd och drar bort i h = min(j+dst_size, height) - j # tar det som är minst j+dst_size eller rasterhöjden och drar bort j # Ovanstående gör att 'sista' raden och kolumnen får en storlek anpassad till vad som är kvar av rasterbilden gdaltranString = "gdal_translate -of GTIFF -srcwin "+str(i)+", "+str(j)+", "+str(w)+", " \ +str(h)+" " + sys.argv[1] + " " + r_dst + "/ruta_"+str(i)+"_"+str(j)+".tif" # skapar en kommandosträng att köra os.system(gdaltranString) # kör kommandosträngen
Några saker om skriptet.
Du behöver ha GDAL installerat! Annars kan du inte köra kommandot ”gdal_translate”. GDAL följer med QGIS men kan installeras separat (http://gdal.org). Prova genom att skriva gdal_translate i ett kommandofönster.
Den första raden är inte Python alls! Det är en kommentar som Python helt struntar i, men i Unix baserade system så talar den om för kommandotolken att skriptet skall exekveras med just Python. I Windows så använder man exempelvis ”python skript.py” för att köra ett Pythonskript och då har man redan i kommandot talat om att det är Python. Om man vill att skript skall kunna köras i flera miljöer (även Mac och Linux/Unix) så bör man dock inkludera denna rad.
Biblioteken os och sys kommer med Python och behövs för att kunna läsa argument och köra kommandon.
Biblioteket gdal behövs för att kunna bearbeta geodatafiler.
Variablerna r_src, dst_size och r_dst tilldelas värden från argument som skickas med skriptet.
Windows python skript.py arg1 arg2 arg3 eller på *nix: ./skript.py arg1 arg2 arg3
Rasterbilden i första argumentet ”mäts” med gdal och bredd och höjd i pixlar lagras i variablerna width och height.
Sedan skapas två loopar i varandra, en för kolumner (i) och en för rader (j). Dessa tilldelar variablerna (i/j) värden från 0 till bildens storlek i pixlar, i steg om det antal pixlar som angivits i det andra argumentet (variabeln dst_size).
Eftersom det kommer att vara sällan en rasterbild går helt jämnt ut uppdelad i den angivna pixelstorleken att dela med, så skapas variablerna w och h, som kompenserar för detta och gör sista kolumnen och raden anpassad i storlek till vad som är kvar.
Det kommando som skall köras i varje loop skapas i textvariabeln gdaltranString.
Här används i, j, w och h med argumentet -srcwin för att plocka ut en ruta i formatet geo-tiff (-of GTIFF) och ger resultatet filnamnet ruta_ kompletterat med rad- och kolumninformation samt filändelsen .tif, lagrat i katalogen r_dst.
Kommandot körs sedan med os.system(gdaltranString).
Om jag tar en godtycklig bild i min ”bilder” katalog på datorn:
Och kör kommandot: ./cut_big_raster.py avatar.jpg 100 avatar. Bilden ”avatar.jpg” är ca 200×200 pixlar stor och det skript jag skapat har jag döpt till cut_big_raster.py.
Resultatfilerna (ruta_0_0.tif, ruta_0_100.tif, ruta_100_0.tif, ruta_100_100.tif) har jag plockat ihop i bilden ovan, men du förstår principen?
Detta skript kunde jag köra på en satellitmosaik över hela Sverige som i original var på 42 Gb! Filen gick att öppna i QGIS, men det tog lite tid. Om jag såg till att vara inzoomad innan jag lade till lagret så var det inga större problem, men med skriptet och en inställning på 5000 pixlar så fick jag i stället 575 st tiffbilder vardera på 75 Mb, vilket i ett virtuellt raster blir mycket lättare att hantera.
I morgon så tänkte jag använda skriptet i QGIS för att skapa ett anpassat verktyg med dialogrutor och liknande för att det skall bli enklare att köra nästa gång jag behöver det.