GDAL/OGR használata Python nyelvben

A nyílt GDAL/OGR (Geospatial Data Abstraction Library, OpenGIS Simple Features Reference Implementation) könyvtárak számos eszközzel segítik a térképészeti adatok feldolgozását. A GDAL a raszteres, az OGR pedig a vektoros adatok kezeléséért felel. A C++ nyelven írt könyvtárakhoz úgynevezett Python kötéseket készítettek, hogy Python programokból elérhetők legyenek a funkciók. A két könyvtár Python kötését a SWIG eszköz segítségével generálják. Az alkalmazásprogramozási felület (API) általában pontosan követi a GDAL és OGR könyvtárak C++ implementációiban használt osztályok és eljárások megnevezéseit. Az elérhető hivatalos dokumentáció (http://gdal.org/python) is automatikusan generált, a meglévő C/C++ függvények leírásai alapján.

Függőségek

A GDAL/OGR könyvtárak használatához a következőkre van szükség:

• Python 2.X (Python 3.X a GDAL 1.7.0 verziójától kezdődően). Python 2-höz ajánlott 2.3 és 2.7 közötti verziót használni.
• libgdal (1.5.0 vagy ennél újabb verzió) a hozzá tartozó header fájlokkal együtt (gdal-devel)
• numpy (1.0.0 vagy ennél újabb verzió) a hozzá tartozó header fájlokkal együtt (numpy-devel)

A numpy telepítése nem feltétlenül szükséges, viszont sok példaprogram és eszköz használja, ezért ajánlott.

Telepítés

GNU/Linux

A GDAL/OGR könyvtár és a Python kötések telepítését elvégezhetjük distutils vagy setuptools segítségével (az utóbbi használata ajánlott). Setuptools-t használva a telepítést a következő paranccsal végezhetjük el:

$ sudo easy_install GDAL

Ha nem szeretnénk feltétlenül a legfrissebb verziót használni és megelégszünk az adott Linux disztribúció repozitóriumában fellelhető verzióval, akkor talán legegyszerűbb a csomagkezelőt használni a telepítésre. Ubuntu operációs rendszeren például a következő paranccsal telepíthető a GDAL/OGR könyvtárak Python kötése:

$ sudo apt-get install gdal-bin python-gdal

Ha forrásból szeretnénk lefordítani a legújabb verziót, akkor a Python kötéshez a --with-python kapcsolót kell használni:

$ ./configure --with-python

Ezután a szokásos módon a make eszköz elvégzi a fordítást és a telepítést:

$ make
$ make install

Windows

A GDAL/OGR könyvtár Python kötésének telepítéséhez a következő lépésekre lesz szükség:

• GDAL Windows bináris fájlok letöltése és kicsomagolása. Alapesetben csak a gdalwin32exe160.zip nevű fájlra lesz szükség. A könyvtárban fellelhető többi fájl bővítmények fejlesztéséhez szükséges. A zip fájlt bárhova ki lehet csomagolni, példa gyanánt a C:\gdalwin32-1.6 könyvtárat fogjuk használni. A kicsomagolás után a PATH rendszerváltozót módosítani kell, hozzá kell adni a C:\gdalwin32-1.6\bin elérési útvonalat.
• Létre kell hozni a GDAL_DATA nevű rendszerváltozót, aminek az értéke az adatokat tartalmazó könyvtár neve, jelen esetben C:\gdalwin32-1.6\data.
• Szükséges lehet az operációs rendszer újraindítása.

Használat

A GDAL Python kötésében öt főmodult érhetünk el a következő módon:

>>> from osgeo import gdal
>>> from osgeo import ogr
>>> from osgeo import osr
>>> from osgeo import gdal_array
>>> from osgeo import gdalconst

Buktatók

Annak ellenére, hogy a Python kötések használata teljesen elrejti a mögöttes C++ könyvtárat, használata pár apró dologban mégis eltér a megszokott Python modulok használatától, amire érdemes figyelnie a Python programozónak:

A Python kötés nem használ kivételeket

A Python kötés alapértelmezetten nem használja a kivételeket, ehelyett meghatározott hibaértéket térít vissza az adott függvény, és egy hibaüzenet jelenik meg a standard kimeneten. Például ha egy nemlétező fájlt próbálunk megnyitni, ez az alapértelmezett működés:

>>> from osgeo import gdal
>>> gdal.Open('/home/user/nincsilyen.img')
ERROR 4: `/home/user/nincsilyen.img' does not exist in the file system,
and is not recognised as a supported dataset name.
>>>

Ha mégis használni szeretnénk a kivételeket, akkor ezeket explicit módon engedélyezni kell az UseExceptions() meghívásával.

>>> from osgeo import gdal
>>> gdal.UseExceptions()        # Kivetelek engedelyezese
>>> gdal.Open('/home/user/nincsilyen.img')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
RuntimeError: `/home/user/nincsilyen.img' does not exist in the file system,
and is not recognised as a supported dataset name.
>>>

Ez megkötés a visszamenőleges kompatibilitás miatt szükséges.

Támogatott adatformátumok

Az OGR könyvtár segítségével számos vektoros formátumot (állománytípust vagy egyéb adatforrást) tudunk kezelni: Aeronav FAA files, AmigoCloud API, ESRI ArcObjects, Arc/Info Binary Coverage, Arc/Info .E00 (ASCII) Coverage, Arc/Info Generate, Atlas BNA, AutoCAD DWG, AutoCAD DXF, CartoDB, Cloudant / CouchDB, CouchDB / GeoCouch, Comma Separated Value (.csv), OGC CSW (Catalog Service for the Web), Czech Cadastral Exchange Data Format, DB2 Spatial, DODS/OPeNDAP, EDIGEO, ElasticSearch, ESRI FileGDB, ESRI Personal GeoDatabase, ESRI ArcSDE, ESRI Shapefile / DBF, FMEObjects Gateway, GeoJSON, Géoconcept Export, Geomedia .mdb, GeoPackage, GeoRSS, Google Fusion Tables, GML, GMT, GPSBabel, GPX, GRASS Vector Format, GPSTrackMaker (.gtm, .gtz), Hydrographic Transfer Format, Idrisi Vector (.VCT), Informix DataBlade, INTERLIS, INGRES, JML, KML, LIBKML, Mapinfo File, Microstation DGN, Access MDB (PGeo and Geomedia capable), Memory, MongoDB, MySQL, NAS - ALKIS, Oracle Spatial, ODBC, MS SQL Spatial, Open Document Spreadsheet, OGDI Vectors (VPF, VMAP, DCW), OpenAir, ESRI FileGDB, OpenStreetMap XML and PBF, PDS, Planet Labs Scenes API, PostgreSQL SQL dump, PostgreSQL/PostGIS, S-57 (ENC), SDTS, SEG-P1 / UKOOA P1/90, SEG-Y, Selafin/Seraphin format, SQLite/SpatiaLite, SUA, SVG, Storage and eXchange Format, UK .NTF, U.S. Census TIGER/Line, VRT - Virtual Datasource, OGC WFS (Web Feature Service), MS Excel format, MS Office Open XML spreadsheet, X-Plane/Flightgear aeronautical data, VDV-451/VDV-452/IDF, Walk, WAsP .map format.

Hozzáférés a filehoz

A különböző fájltípusoknak és más adatforrásoknak a kezelésére az OGR könyvtár úgynevezett meghajtókat (vagy driver­eket) használ. A következő Python kóddal tudjuk megvizsgálni, hogy milyen driver­ek állnak rendelkezésünkre:

from osgeo import ogr
driverList = []

for i in range(ogr.GetDriverCount()):
  driver = ogr.GetDriver(i)
  driverName = driver.GetName()
  if not driverName in driverList:
    formatsList.append(driverName)

for i in formatsList:
  print i

A kiírt nevek alapján azonosíthatjuk a megfelelő drivert. A Shape fájl formátumot például az „ESRI Shapefile” nevű meghajtóval kezelhetjük. A megfelelő driver név szerint is elérhető (ha nincs telepíve az adott nevű driver, akkor a GetDriverByName függvény None értéket térít vissza). A Shape fájlokat kezelő meghajtót tehát a következő függvényhívással érhetjük el: driver = ogr.GetDriverByName('ESRI Shapefile') Az állományt ezután a meghatón keresztül nyitjuk meg az Open függvény segítségével, aminek első paramétere az állomány neve (teljes elérési útvonal), a második pedig egy egész szám, aminek értéke 0 vagy 1 lehet. A második paraméter 0, ha az állományt csak olvasásra nyitjuk meg, az érték 1 ha írni is szeretnénk bele.

file = driver.Open(filename, 0)
if file is None:
  print ('Nem tudtam megnyitni a fájlt!')

Amennyiben a meghajtó nem tudta megnyitni az állományt, akkor a None értéket adja vissza. Ez a helyzet akkor fordulhat elő, ha a Shape fájl tartalma sérült vagy az shx vagy dbf fájl nem található.

Térképészeti adat kinyerése

A következő lépés a Shape fájlban található réteghez (layer) valő hozzáférés. Ezt a funkciót a

GetLayer(index) függvény biztosítja. Shape fájlok esetében az index mindig 0 (vagy el is lehet hagyni ezt a paramétert), az index csak olyan formátumok esetében hasznos, mint pl. a GML vagy a TIGER. A következő sorral tehát a Shape fájl egyetlen rétegét szerezzük be:

layer = datasource.GetLayer()

Ezután következik a rétegen található elemek (features) beolvasása. A feature­ek számát a layer GetFeatureCount() függvényével kérhetjük le, és az egyes feature­ek a GetFeature(index) függvénnyel érhetjük el. Vagy végig lehet menni az összes feature­en a következő kóddal:

feature = layer.GetNextFeature()
while feature:
  # feldolgozás
  feature = layer.GetNextFeature()
  layer.ResetReading() #ha újra kell kezdeni a beolvasást

A Shape fájlunk csak egyetlen feature­t tartalmaz, ezért a GetNextFeauter() egyszeri meghívásával megoldjuk a hozzáférést. Az elem mértani objetkumát a GetGeometryRef() függvénnyel kérhetjük le, típusát pedig a GetGeometryType() vagy GetGeometryName() függvénnyekkel ellenőrizhetjük le. A következő sorokban például azt ellenőrizzük le, hogy a poligon vagy multipoligon típusú elemmel van­e dolgunk:

geometry = feature.GetGeometryRef()
if geometry.GetGeometryName() == 'POLYGON' or geom.GetGeometryName() =='MULTIPOLYGON':
  # feldolgozás

A típusokat az OGR konstansaival (pl. ogr.wkbPoint, ogr.wkbLineString, ogr.wkbPolygon, ogr.wkbMultiPoint, ogr.wkbMultiLineString, ogr.wkbMultiPolygon, stb.) is azonosíthatjuk.

Attribútumok

Az attribútomokat a GetField() függvénnyel és annak variációival érhetjük el. Például:

attr = feature.GetField('id')
attrstr = feature.GetFieldAsString('id')

Az attribútumok számát a GetFieldCount() függvénnyel kapjuk meg.

Példaprogramok

Shape file geometriájának egyszerűsítése

A következő Python kód egy Shape file geometriáját egyszerűsíti, adott toleranciával. A fontosabb sorokhoz magyarázatot fűztünk.

#!/usr/bin/python
# ­*­ coding: utf­8 ­*­

import os, sys
from osgeo import ogr

# infile ­ bemeneti állomány neve
# outfile ­ kimeneti állomány neve
# tolerance ­ a egyszerűsítés paramétere

def simplify(infile, outfile, tolerance):
  # az ESRI Shapefile meghajtó
  driver = ogr.GetDriverByName('ESRI Shapefile')
  
  # olvasásra nyitjuk meg a bemeneti állományt
  infile = driver.Open(infile,0)
  if infile is None:
    print 'Nem tudom megnyitni a(z) ', infile, ' nevű állományt!'
    sys.exit(1)

  # a bemeneti állomány adatai: layer, feature, geometry
  inputLayer = infile.GetLayer()
  inputFeature = inputLayer.GetNextFeature()
  geom = inputFeature.GetGeometryRef()
  geomType = geom.GetGeometryType()

  # a kimeneti állomány létrehozása
  if os.path.exists(outfile):
    os.remove(outfile)
  try:
    output = driver.CreateDataSource(outfile)
  except:
    print 'Nem tudom létrehozni a(z)', outfile, ' nevű állományt!'
    sys.exit(1)

  # réteg létrehozása a kimeneten
  outputLayer = output.CreateLayer('Tolerance',geom_type=geomType,srs=inputLayer.GetSpatialRef())
  if outputLayer is None:
    print 'Nem tudom lérehozni a megfelelő layer­t a kimeneti állományban!'
    sys.exit(1)
  outputLayerDef = outputLayer.GetLayerDefn()
  featureID = 0

# végigmegyünk az összes elemen
while inputFeature:
  # az eredeti geometria
  geometry = inputFeature.GetGeometryRef()
  # az egyszerűsített geometria
  simplifiedGeom = geometry.Simplify(tolerance)
  # megpróbálunk létrehozni egy új feature­t az egyszerűsített geometriával
  try:
    newFeature = ogr.Feature(outputLayerDef)
    newFeature.SetGeometry(simplifiedGeom)
    newFeature.SetFID(featureID)
    outputLayer.CreateFeature(newFeature)

  except:
    print "Nem tudtam létrehozni az egyszerűsített geomatriát!"

  newFeature.Destroy()
  inputFeature.Destroy()
  inputFeature = inputLayer.GetNextFeature()
  featureID += 1
  infile.Destroy()
  output.Destroy()
  print "Az egyszerűsítést sikeresen végrehajtottam"

return

# a parancssorban átadott paraméterek
if (len(sys.argv) < 4):
  print ('Használat: python simplify.py <bemenet> <kimenet> <tolerancia>')
  sys.exit(0)

# az egyszerűsítő függvény meghívása
simplify(sys.argv[1], sys.argv[2], float(sys.argv[3]))

Népsűrűség ábrázolása térképen

Az alábbi kód, egy meglévő shapefile rétegre egy másikat generál, egy megadott atríbútum táblából, aminek segítségével szemlélteti, hogy az adott területeken mekkora a népsűrűség. Ehhez pontokat generál, egy pont 100 embert reprezentál.

from osgeo import ogr
import random
# shapefile megnyitasa ogr reteg letrehozasa elso feature lekerese
source = ogr.Open("GIS_CensusTract_poly.shp")
county = source.GetLayer("GIS_CensusTract_poly")
feature = county.GetNextFeature()
# kimeneti shapefile es reteg letrehozasa
driver = ogr.GetDriverByName('ESRI Shapefile')
output = driver.CreateDataSource("PopDensity.shp")
dots = output.CreateLayer("PopDensity", geom_type=ogr.wkbPoint)
while feature is not None:
  field_index = feature.GetFieldIndex("POPULAT11")
  population = int(feature.GetField(field_index))
  # 1 pont szaz embert reprezental 
  density = population / 100
  count = 0   
  while count < density:
    geometry = feature.GetGeometryRef()
    minx, maxx, miny, maxy = geometry.GetEnvelope()
    x = random.uniform(minx,maxx)
    y = random.uniform(miny,maxy)
    f = ogr.Feature(feature_def=dots.GetLayerDefn())
    wkt = "POINT(%f %f)" % (x,y)
    point = ogr.CreateGeometryFromWkt(wkt)
    # Csak akkor hasznaljuk a pontot, ha benne van az adott alakzatban
    if feature.GetGeometryRef().Contains(point):
        f.SetGeometryDirectly(point)
        dots.CreateFeature(f)
        count += 1
    # objektum eltorlese
    f.Destroy()
  feature = county.GetNextFeature()
source.Destroy()
output.Destroy()

Hivatkozások

http://www.gdal.org/

https://pcjericks.github.io/py-gdalogr-cookbook/

http://www.osgeo.org/gdal_ogr

http://geoexamples.com/

http://www.digital-geography.com/create-and-edit-shapefiles-with-python-only/