Tag Archive | gdal

Skapa geodata från Python

Igår använde jag ExifRead för att läsa GPS taggar från foton, men för att skapa punkter baserade på dessa data så krävs det något mera. I ArcMap så finns det inbyggda funktioner för detta, men nu tänkte jag se om det går att göra detta helt i Python.

Läs mer…

Annonser

GDAL

Som framgått i sociala media för ett tag sedan så har man uppdaterat QGIS till GDAL 2.0 (från 1.11). Vad betyder det för dig då?

Jo exempelvis så öppnar det för möjligheten att läsa och redigera MapInfo TAB filer i QGIS. (Även MIF skall fungera.)

Om detta är en stor sak eller ej vet jag inte, då jag inte använder MapInfo och sällan stöter på TAB-filer, men det finns fler fördelar med uppdatering av GDAL i QGIS. Läs mer…

ArcToolbox för GDAL vrt

Stor glädje! Jag har fått min toolbox för att skapa GDAL vrt-filer i ArcMap att fungera!

Använder du ArcMap och vill prova så kan du bläddra ner mot slutet av artikeln, där finns en länk till de filer du behöver ladda hem. Läs mer…

Skapa GDAL vrt

Jag följer upp gårdagens artikel om virtuella rasterkataloger i ArcMap med en metod för att skapa denna typ av kataloger från terminalen/kommandoprompten.

Jag antar att man sitter på en Windowsmaskin när man använder ArcMap så därför är instruktionen anpassad för Windows. När du läst klart denna korta artikel så har du allt du behöver för att skapa vrt-filer även för ArcMap, och allt du behöver är GDAL. Läs mer…

GDAL och Python

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

Skärmbild från 2015-03-07 10:11:08Det 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:

Skärmbild från 2015-03-07 10:47:01GDAL 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…

Förenkla med skript

Många arbetsuppgifter görs flera gånger och inte nödvändigtvis så ofta att man alltid kommer ihåg hur man gjorde. Då kan det vara läge att fundera över att skapa ett skript som löser uppgiften åt dig.

Jag använder ibland GPS när jag är ute och gör något, och det är då smidigt att ha ett skript som automatiskt tankar av och tömmer mottagaren till datorn, samt gör om data till ESRI shape. Jag har visat detta tidigare i bloggen, men jag tänkte vara lite mera generell denna gången.

Har du ett ofta återkommande behov av att flytta, spara, konvertera eller bearbeta geografiska data så är detta inlägget för dig.

Jag kommer att bygga mitt exempel på terminalen i Linux och ett verktyg i biblioteket GDAL, som även finns för Mac och Windows.

Vad är det som är så bra med terminalen? Jo alla kommandon som går att skriva i terminalen kan enkelt skriptas på ett enkelt och interaktivt sätt utan att man behöver avancerad programmering eller ens tillverka ett grafiskt användargränssnitt, vilket inte alltid är så enkelt. Du kan exempelvis enkelt skapa ett skript som söker efter och kopierar alla QGIS projektfiler till en gemensam katalog eller skapar en zip-fil med alla Jpg-bilder på en bestämd sökväg.

För att hantera geografiska data i olika format så använder jag verktyget ogr2ogr som installeras med GDAL via kommandot:

sudo apt-get install gdal-bin

Vill du se vad kommandot kan utföra så skriv ogr2ogr –help och vill du se alla format som verktyget kan hantera så skriv ogr2ogr –formats.

”ESRI Shapefile” (read/write)
”MapInfo File” (read/write)
”UK .NTF” (readonly)
”SDTS” (readonly)
”TIGER” (read/write)
”S57” (read/write)
”DGN” (read/write)
”VRT” (readonly)
”REC” (readonly)
”Memory” (read/write)
”BNA” (read/write)
”CSV” (read/write)
”GML” (read/write)
”GPX” (read/write)
”KML” (read/write)
”GeoJSON” (read/write)
”Interlis 1” (read/write)
”Interlis 2” (read/write)
”GMT” (read/write)
”SQLite” (read/write)
”DODS” (readonly)
”ODBC” (read/write)
”PGeo” (readonly)
”OGDI” (readonly)
”PostgreSQL” (read/write)
”MySQL” (read/write)
”PCIDSK” (readonly)
”XPlane” (readonly)
”AVCBin” (readonly)
”AVCE00” (readonly)
”DXF” (read/write)
”Geoconcept” (read/write)
”GeoRSS” (read/write)
”GPSTrackMaker” (read/write)
”VFK” (readonly)

Jag tänker göra ett skript som tar gpxfiler och adderar dessa data till två shapefiler. En för spårloggar och en för waypoints. Det skulle lika gärna kunnat vara en PostGIS databas eller en GeoJSON fil.

Att bara göra om filer, i det här fallet från gpx till shp, är inte alls svårt. Använd bara kommandot ogr2ogr -f ”ESRI Shapefile” målfil.shp källfil.gpx ”lager”. I praktiken så blir det dock lite mera, inte minst då jag vill skriva till samma fil varje gång.

Nåja låt oss sätta igång. Börja med att skapa en textfil och döp den till ”gps2shp.sh” eller valfritt filnamn med filändelsen *.sh. Först i filen skriver man #!/bin/bash för att tala om att det är ett skript som använder ”bash”.

Sedan vill jag att användaren skall mata in eller dra-och-släppa gpxfilen, så då skriver man följande för att hantera detta:

echo -n ”Skriv in sökväg till (eller dra-och-släpp) gpx-fil >”
read Filnamn
Filnamn=${Filnamn//\’/}

Efter en bekräftelse med retur så finns nu filens sökväg lagrad i variabeln ”Filnamn”. Den sista raden tar bort extratecken som skapas vid ”dra-och-släpp”, annars fungerar det inte.

Man bör även sätta projektionen på shapefilen så detta kan man skapa en variabel för och använda senare:

Prj=”-a_srs EPSG:4326”

Konverteringen från gpx till shape får även datum och tid att krångla så ett kommando för att konvertera datum till text kan också sparas i en variabel:

DatumFix=”-fieldTypeToString DateTime”

För att hålla isär saker och ting så kan man även skapa en variabel för sökvägen till målfilerna:

Target=”gps-loggar”

I och med att det blir olika kommandon om loggen skall skapas eller om den redan finns så måste man använda ett villkor i skriptet. Villkoret testar om det finns en känd målfil på angiven sökväg.

Testa=$Target”/waypoints.shp”
if [ -f $Testa ]
then
#Kod att köra om filen finns
else
#Kod att köra om filen inte finns
fi

Jag återkommer till den utelämnade koden längre ner. Kommandot -f $Testa undersöker om filen i variabeln ”Testa” existerar, och om den gör det är villkoret sant. Om filen inte finns är villkoret falskt.

Det kommando som skall köras om filen finns ryms på en enda rad:

ogr2ogr -update -append $Target $Filnamn ”tracks” ”waypoints” $DatumFix $Prj

Koden som behöver köras om filen inte finns är snarlik, men föregås av att sökvägen till loggen skapas:

mkdir $Target
ogr2ogr $Target $Filnamn ”tracks” ”waypoints” $DatumFix $Prj

i Båda kommandona ovan så väljs endast lagerna tracks (linjer) och waypoints (punkter) ut. Vill du se vilka fler lager som finns i en fil skriver du in kommandot ogrinfo filnamn.

1: waypoints (Point)
2: routes (Line String)
3: tracks (Multi Line String)
4: route_points (Point)
5: track_points (Point)

Det var det hela. Man kan ju förstås lägga till en rad för att tala om att allt skriptet är klart.

clear
echo ”Filen $Filnamn lagrades på sökvägen $Target”
sleep 5

För att kunna köra skriptet när man dubbelklickar på filen så kryssar man i ”Tillåt körning av filen som ett program” på fliken ”Rättigheter” under filegenskaperna (högerklicka på filen).

Skärmbild från 2013-08-24 19:38:57

Nu är det bara att köra skriptet när man behöver lägga till gpx filer till en loggfil, och kom ihåg att principen är den samma oavsett vilka format du vill konvertera emellan så länge de stöds av kommandot. Loggfilerna skapas där skriptet är sparat som skriptet är skrivet nu, men det går att ange en annan sökväg enkelt genom att redigera skriptet. Kanske även lägga till en ”dra-och-släpp” funktion för att ange sökvägen till var loggarna skall sparas?

Du kan även hämta hela skriptet på sökvägen http://geosupport.kvarnarp.eksjo.com/files/gpx2shp.sh.