I detta inlägg tänkte jag introducera GDAL och hantering av geodata i Python.
GDAL (http://gdal.org) är ett översättningsbibliotek som gör det möjligt att hantera såväl vektor- som rasterdata i flera olika programmeringsspråk däribland Python.
Detta öppnar flera intressanta möjligheter för att hantera data i olika sammanhang. Om du har en metod att skriva Python skript så kan du med biblioteken hantera geodata, och det behöver inte vara exempelvis ArcGIS eller QGIS (som båda hanterar Python), utan som jag skall visa här så räcker det med ett terminalfönster.
Först så gäller det att ha GDAL/OGR installerat (GDAL hanterar rasterdata och OGR vektordata). Hur du gör detta beror på vilket operativsystem du har så jag hoppar över det här. Gå till gdal.org och läs mer där.
Chansen är att du redan har det installerat så börja med att öppna en terminal eller ett kommandofönster och skriv in:
python from osgeo import ogr
Det första kommandot startar en Python-tolk och det andra läser in ett bibliotek som hanterar vektordata.
Så länge du inte stöter på något felmeddelande så skall det bara vara att köra vidare.
För att kunna hantera sökvägar så behövs ännu ett bibliotek:
import os
Jag tänkte ta reda på hur många grillplatser det finns i Halmstad (ett gammalt lager jag råkar ha liggande) så jag börjar med att skapa en variabel som innehåller sökvägen till shapefilen:
grillplatser = r"/home/klakar/geodata/Grillplatser.shp"
Du får naturligtvis anpassa detta kommando till en shapefil du har tillgång till, men det kan i princip vara vilken typ av lager som helst. Bokstaven ”r” i början anger att texten är skriven i ”raw” format. Det som står innanför citationstecken skall med andra ord endast tolkas som text och inte på något annat sätt. Jämför nedanstående kommandon:
print "Rad ett\nRad två" print r"Rad ett\nRad två"
Sökvägen skall läsas in som datakälla och då behöver man även tala om vilket format som skall användas för att tolka dessa data:
drivrutin = ogr.GetDriverByName('ESRI Shapefile')
dataKalla = drivrutin.Open(grillplatser, 0)
Den sista nollan anger att filen skall läsas men inte gå att skriva till. En etta gör filen skrivbar.
God sed är att testa att filen finns och gick att läsa in:
if dataKalla is None: print 'Kan inte öppna %s' % (grillplatser)
Här ovan så händer lite intressanta saker för den som inte är van vid Python. Det första kommandot anger att ”nu kommer det mera” och därför exekveras inget vid varje ny rad. Detta indikeras av de tre punkterna (…) i början på varje rad. Om det inte finns något i variabeln dataKalla så skall det efterföljande köras…
En annan sak som är viktig att känna till om Python är att indrag är viktiga. För att tala om för Python vilken kod som skall köras om vilkoret i if-satsen är sant så skall det finnas ett indrag vid alla dessa rader. Den första efterföljande raden som inte har indrag, är den första kodrad som körs om vilkoret är falskt:
else: print 'Öppnade %s' % (grillplatser) lager = dataKalla.GetLayer() antalPlatser = lager.GetFeatureCount() print "Antal grillplatser i %s: %d" % (os.path.basename(grillplatser),antalPlatser)
Om det gick bra att läsa in lagret så körs ”else” i stället för ”if”.
Här, och i förra kodraderna, så står det %s med mera. Detta är ett sätt i Python att skapa en platshållare för något som definieras efter raden. %s ersätts i ett fall av innehållet i variabeln grillplatser.
Lager är en variabel där man läser in datalagret från den definierade datakällan och sedan är det bara att använda ogr-funktionen för att räkna antalet objekt (GetFeatureCount) och lägga in resultatet i en ny variable.
Den sista raden skriver ut resultatet till skärmen och här används två %-platshållare.
Hela skriptet:
GDAL och OGR kan nu vansinnigt mycket mera! Det går att öppna ett vektorlager, skapa en buffer och spara som ett nytt lager, bara med några rader kod och det utan att över huvud taget ha öppnat ett GIS-program.
Det går att skapa data från grunden, konvertera format (137 rasterformat och 81 vektor), med mera.
Det finns terminalkommandon som är byggda för att förenkla några av dessa funktioner.
ogr2ogr – konvertera mellan olika vektorformat.
pgsql2shp – shp2pgsql – konvertera mellan postgresql och shape.
Har du egna format du ofta konverterar mellan så är det inte så svårt att skapa ett liknande kommando själv, med hjälp av GDAL och OGR. Det kräver dock lite kunskaper i programmering, men det kanske är dags att lära sig lite Python…