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. 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.

Telepítés

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.

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

A két könyvtár Python kötését a SWIG eszköz segítségével generálják. 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.

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.

A Python összeomlik, ha egy olyan objektumot használunk, ami összefügg egy már előzőleg törölt objektummal

A GDAL és OGR objektumok valójában C++ objektumok, amik mutatókkal kapcsolódnak egymáshoz. Megtörténhet, hogy a Python garbage collector felszabadítja egy objektum tárhelyét (vagy a programozó töröl egy objektumot del-lel), miközben egy másik objektum mutatóval hivatkozik rá. Ilyenkor a Python összeomlik. A következő példában egy bemeneti állományt megnyitva létrehozunk egy objektumot, majd egy másikat, ami a háttérben hivatkozik az elsőre (C++ mutatóval). Ezután kitöröljük az első objektumot és meghívjuk a másodiknak egy metódusát. Ennyi elég ahhoz, hogy a Python összeomoljon.

>>> from osgeo import gdal
>>> adat = gdal.Open('/home/user/valami.img')
>>> band = adat.GetRasterBand(1) # Egy uj objektum letrehozasa, ami kapcsolatban van az adat objektummal
>>> del adat           # Toroljuk az adat objektumot
>>> band.GetMinimum();    # A band objektum egy metodusanak hivasa es a Python osszeomlik
Segmentation fault (core dumped)

A következő példában a garbage collector törli automatikusan az osgeo.gdal.Dataset típusú objektumot még mielőtt a GetMinimum() metódus meghívódna, mert a osgeo.gdal.Band objektum létrehozása után erre már nincs szükség (a Python "logikája" szerint, ami nem tudhat a háttérben levő C++ implementáció részleteiről).

>>> from osgeo import gdal
>>> print gdal.Open('/home/user/valami.img').GetRasterBand(1).GetMinimum()
Segmentation fault (core dumped)

Fájlkezelés

A különböző fájltípusoknak és más adatforrásoknak a kezelésére az GDAL/OGR könyvtárak úgynevezett meghajtókat (vagy driver­eket) használnak. Például a következő Python kóddal tudjuk megvizsgálni, hogy milyen OGR 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ámogatott adatformátumok

A GDAL segítségével számos raszteres formátumot (állománytípust vagy egyéb adatforrást) tudunk kezelni. Az OGR könyvtár vektoros formátumú fájlokat kezel. A következőkben felsoroljuk a támogatott raszteres illetve vektoros állománytípusokat:

• Raszteres adatformátumok (GDAL): Arc/Info ASCII Grid, ACE2, ADRG/ARC Digitilized Raster Graphics (.gen/.thf), Arc/Info Binary Grid (.adf), AIRSAR Polarimetric, Azavea Raster Grid, Magellan BLX Topo (.blx, .xlb), Bathymetry Attributed Grid (.bag), Microsoft Windows Device Independent Bitmap (.bmp), BSB Nautical Chart Format (.kap), VTP Binary Terrain Format (.bt), CALS Type I, CEOS (Spot for instance), TerraSAR-X Complex SAR Data Product, USGS LULC Composite Theme Grid, DirectDraw Surface, Spot DIMAP (metadata.dim), DODS / OPeNDAP, First Generation USGS DOQ (.doq), New Labelled USGS DOQ (.doq), Military Elevation Data (.dt0, .dt1, .dt2), Arc/Info Export E00 GRID, ECRG Table Of Contents (TOC.xml), ERDAS Compressed Wavelets (.ecw), ESRI .hdr Labelled, Erdas Imagine Raw, NASA ELAS, ENVI .hdr Labelled Raster, Epsilon - Wavelet compressed images, ERMapper (.ers), Envisat Image Product (.n1), EOSAT FAST Format, FITS (.fits), Generic Binary (.hdr Labelled), Oracle Spatial GeoRaster, GSat File Format, Graphics Interchange Format (.gif), WMO GRIB1/GRIB2 (.grb), GMT Compatible netCDF, GRASS Raster Format, GRASS ASCII Grid, Golden Software ASCII Grid, Golden Software Binary Grid, Golden Software Surfer 7 Binary Grid, Generic Tagged Arrays (.gta), TIFF / BigTIFF / GeoTIFF (.tif), GXF - Grid eXchange File, Hierarchical Data Format Release 4 (HDF4), Hierarchical Data Format Release 5 (HDF5), HF2/HFZ heightfield raster, Erdas Imagine (.img), Image Display and Analysis (WinDisp), ILWIS Raster Map (.mpr,.mpl), Intergraph Raster, IRIS, ISCE raster, USGS Astrogeology ISIS cube (Version 2), USGS Astrogeology ISIS cube (Version 3), JAXA PALSAR Product Reader (Level 1.1/1.5), Japanese DEM (.mem), JPEG JFIF (.jpg), JPEG-LS, JPEG2000 (.jp2, .j2k), JPEG2000 (.jp2, .j2k), JPEG2000 (.jp2, .j2k), JPEG2000 (.jp2, .j2k), JPEG2000 (.jp2, .j2k), JPIP (based on Kakadu), KEA, KRO, NOAA Polar Orbiter Level 1b Data Set (AVHRR), Erdas 7.x .LAN and .GIS, FARSITE v.4 LCP Format, Daylon Leveller Heightfield, MBTiles, OziExplorer .MAP, In Memory Raster, Vexcel MFF, Vexcel MFF2, MG4 Encoded Lidar, Meta Raster Format, Multi-resolution Seamless Image Database, Meteosat Second Generation, EUMETSAT Archive native (.nat), NLAPS Data Format, NOAA NGS Geoid Height Grids, NITF (.ntf, .nsf, .gn?, .hr?, .ja?, .jg?, .jn?, .lf?, .on?, .tl?, .tp?, etc.), NetCDF, OGDI Bridge, OZI OZF2/OZFX3, PCI .aux Labelled, PCI Geomatics Database File, PCRaster, Geospatial PDF, NASA Planetary Data System, Planet Labs Mosaics API, Portable Network Graphics (.png), Netpbm (.ppm,.pgm), R Object Data Store, Rasdaman, Rasterlite - Rasters in SQLite DB, Swedish Grid RIK (.rik), Raster Matrix Format (*.rsw, .mtw), ROI_PAC Raster, Raster Product Format/RPF (CADRG, CIB), RadarSat2 XML (product.xml), Idrisi Raster, Sentinel 1 SAR SAFE (manifest.safe), Sentinel 2, SAGA GIS Binary format, SAR CEOS, ArcSDE Raster, USGS SDTS DEM (*CATD.DDF), SGI Image Format, Snow Data Assimilation System, Standard Raster Product (ASRP/USRP), SRTM HGT Format, Terragen Heightfield (.ter), USGS ASCII DEM / CDED (.dem), OGC Web Coverage Service, WEBP, OGC Web Map Service, and TMS, WorldWind, On Earth tiled, VirtualEarth, ArcGIS REST, IIP), OGC Web Map Tile Service, X11 Pixmap (.xpm), ASCII Gridded XYZ, ZMap Plus Grid, here

• Vektorgrafikus adatformátumok (OGR): 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.

Példaprogramok

Térképészeti adat kinyerése OGR-ban

Röviden megmutatjuk az ESRI 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.

OGR 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.

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/