Εισαγωγή

H GDAL (Geospatial Data Abstraction Library) είναι μια ελεύθερη βιβλιοθήκη (free software),που διαχειρίζεται γεωγραφικά δεδομένα. Η GDAL αποτελείται από δυο τμήματα: την GDAL που διαχειρίζεται raster δεδομένα και την OGR, η οποία διαχειρίζεται vector. Λόγω της ευρείας χρήσης της έχει μπει κάτω από την ομπρέλα του OSGeo (Open Source Geospatial Foundation). Η δύναμη της είναι τέτοια, που χρησιμοποιείται σε πολλά πακέτα GIS, τόσο ελέυθερου όσο και κλειστού/εμπορικού λογισμικού. Μερικά από αυτά είναι: Grass GIS, Quantum GIS, Google Earth, MapServer και ArcGIS.

Η βιβλιοθήκη είναι γραμμένη σε C, αλλά έχει μεταφερθεί και σε Python.

Εγκατάσταση σε Ubuntu 10.04

Καλό θα είναι να εγκατασταθεί όλο το αποθετήριο του UbuntuGis που παρέχει άμεση πρόσβαση σε διάφορα προγράμματα γεωγραφικού ενδιαφέροντος. Για να γίνει αυτό ανοίξτε ένα terminal (Applications → Accessories → Terminal) και δώστε τις παρακάτω εντολές:

  sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
  sudo apt-get update

Από δω και στο εξείς μπορείτε να προσθέτετε πακέτα από το UbuntuGis κατευθείαν από τον Synaptic Package Manager. Για να προσθέσετε την GDAL ανίξτε το Synaptic Package Manager και ψάξτε το αρχείο gdal-bin. Επίσης για χρήση στην Python προσθέστε και το αρχείο python-gdal.

GDAL module

Όπως αναφέρθηκε, το module διαχειρίζεται raster αρχεία, δηλαδή με απλά λόγια είναι χρήσιμο για την Ψηφιακή Τηλεπισκόπηση! Τα είδη των αρχείων που υποστηρίζονται περιλαμβάνουν όλα τα γνωστά αρχεία εικόνα (TIFF, Erdas, ERmapper κ.λπ.) και πολλά άλλα! Μια πλήρη λίστα με αυτά υπάρχει εδώ!

Module στην Python δεν είναι τίποτα άλλο από ένα αρχείο .py (αρχείο Python δλδ), το οποίοι περιέχει μέσα έτοιμες συναρτήσεις και κλάσεις, οι οποίες καλώντας τες κάνουν χρησιμα πράγματα. Τα modules πρέπει να κλήθούν στην αρχή του κάθε προγράμματος ώστε να είναι διαθέσιμο το περιεχόμενό τους. Για παράδειγμα, το math είναι ένα μαθηματικό module της Python και έχει μέσα διάφορες συναρτήσεις, οι οποίες είναι έτοιμες για χρήση και μας γλιτώνουν από υπολογισμούς (πχ η μετατροπή από rad σε degrees γίνεται από την συνάρτηση degrees(x)).

Το GDAL module αρθρώνεται σε κάποιες αυτόνομες συναρτήσεις και κάποιες κλάσεις, οι οποίες έχουν μέσα τους άλλες συναρτήσεις.

Άνοιγμα αρχείου

Προχωράμε αμέσως στο άνοιγμα ενός αρχείου. Έστω ότι έχουμε μια εικόνα .ers (ERMapper format) στο ‘/home/user/Desktop’ με ονομα ‘image1.ers’. Αυτό γίνεται σε 3 στάδια:

⇒ Κλήση των modules της GDAL.

⇒ Ενεργοποίηση των οδηγών (drivers) που είναι κατάλληλοι για το άνοιγμα του αρχείου .ers

⇒ Άνοιγμα αρχείου

Έτσι έχουμε:

  from osgeo import gdal
  from osgeo.gdalconst import *
  gdal.AllRegister()
  image=gdal.Open('/home/user/Desktop/image1.ers', GA_ReadOnly)

Στο συγκεκριμένο παράδειγμα επειδή δεν θέλουμε να επεξαργαστούμε την εικόνα αλλά μόνο να διαβάσουμε τα δεδομένα της χρησιμοποιήσαμε έναν driver-’μπαλαντέρ’ (τον AllRegister()), ο οποίος ανοίγει όλες τις εικόνες, μόνο όμως για ανάγνωση.

Ορίζουμε στην μεταβλητή image την διαδικασία ανοίγματος της εικόνας. Η εικόνα ανοίγει χρησιμοποιώντας την μέθοδο Open() της κλάσσης gdal, η οποία δέχεται 2 arguments: το path της εικόνας και την κατάσταση ανοίγματος.

Προσπέλαση Βασικών Δεδομένων

Μπορούμε να πάρουμε κάποιες γενικές πληροφορίες για το αρχείο ως εξής:

Η συνάρτηση GetGeoTransform() μας επιστρέφει ένα tuple με τις συν/νες του πάνω αριστερά pixel, την γωνία στροφής της εικόνας κατα Χ και Υ άξονα καθώς και το μέγεθος του pixel (πάλι) κατά Χ και Υ άξονα. To tuple είναι ΠΑΝΤΑ της μορφής (Χ, μέγεθος pixel κατα Χ, στροφή κατα Χ, Υ, στροφή κατά Υ, Μέγεθος pixel κατά Υ). Παρατίθεται ένα πραγματικό output για μια τυχαία εικόνα:

  (219543.05899931892, 30.0, 0.0, 4407229.9912246773, 0.0, -30.0)

Η συνάρτηση GetTProjection() μας επιστρέφει ένα string με τα στοιχεία προβολής της εικόνας. Το output είναι αυτής της μορφής (για μια τυχαία εικόνα):

  'PROJCS["EGSA87",PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",23.99999882666041],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

Τέλος οι μεθόδοι RasterXSize, RasterYSize και RasterCount μας επιστρέφουν τον αριθμό των pixels στον Χ και Υ άξονα και τον αιρθμό των καναλιών της εικόνας.

Από τα παραπάνω μπορούμε να συνθέσουμε ένα απλό πρόγραμμα, το οποίο ανοίγει μια εικόνα και μας επιστρέφει τα βασικά της στοιχεία:

  from osgeo import gdal
  from osgeo.gdalconst import *
  gdal.AllRegister()
  image=gdal.Open('/home/user/Desktop/image1.ers', GA_ReadOnly)
  print 'Pixel ston X aksona:', image.RasterXSize
  print 'Pixel ston Y aksona:', image.RasterYSize
  print 'Arithmos Kanaliwn:', image.RasterCount 
  print 'Geodetiko Datum:', image.GetProjection().split('"')[1]
  print 'Provoli:', image.GetProjection().split('"')[3]
  print 'Syntetagmenes panw aristerou pixel:', image.GetGeoTransform()[0], ',', image.GetGeoTransform()[3]
  print 'Strofi kata X kai kata Y:', image.GetGeoTransform()[2],',',image.GetGeoTransform()[4]
  print 'Diastaseis pixel kata X kai Y:', image.GetGeoTransform()[1],',',math.fabs(image.GetGeoTransform()[5])

OGR module

To OGR είναι το αντίστοιχο module διαχείρισης vector δεδομένων.

Άνοιγμα αρχείου

Το άνοιγμα αρχείων γίνεται σχεδόν όπως στην GDAL. Για παραδείγμα εδώ θα χρησιμοποιήσουμε ένα αρχείο Shapefile (.shp). Ο λόγος είναι ότι το shapefile είναι ένα από τα πιο διαδεδομένα αρχεία vector δεδομένων. Έστω λοιπόν ότι έχουμε ένα αρχείο test.shp. Αφήνω προς το παρόν την γεωμετρία του (point, line, polygon) φλου, επειδή θα αναφερθούμε σε εντολές που ισχύουν για όλα τα παραπάνω είδη. Έτσι έχουμε:

  from osgeo import ogr
  driver=ogr.GetDriverByName('ESRI Shapefile')
  datasource=driver.Open('/home/user/Desktop/test.shp', 0)
  layer=datasource.GetLayer()

⇒ Line 1 : καλώ το ogr module

⇒ Line 2 : στην περίπτωση του ogr πρέπει να καλέσω συγκεκριμένο driver. O driver που καλώ είναι ο ‘ESRI Shapefile’ (προσοχή στα κεφαλαία/μικρά). Τον καλώ μέσω της μεθόδου GetDriverByName και τον αποθηκεύω στην μεταβλητή driver.

⇒ Line 3 : ανοίγω το αρχείο όπως και στην περίπτωση της GDAL. Η μέθοδος Open() δέχεται εκτός από το path του αρχείου και μια τιμή, η οποία έιναι 0 ή 1. Το 0 χρησιμοποιείται για read-only άνοιγμα, ενώ το 1 σε περίπτωση που θέλουμε να ‘πειράξουμε’ το αρχείο.

⇒ Line 4 : Η βασική πληροφορία του Shapefile βρίσκεται στο Layer. Συνήθως οι βασικές εργασίες GIS γίνονται στο επίπεδο Layer ή στο επίπεδο Feature.

Πειγραφή του Shapefile

Το shapefile είναι ένα διανυσματικό αρχείο γεωγραφικής πληροφορίας. Το κάθε αρχείο απεικονίζει αποκλειστικά ένα είδος πληροφορίας (σημειακή, γραμμική ή επιφανειακή). Παρόλο που συχνά το shapefile παρουσιάζεται ως αρχείο, στην πραγματικότητα είναι μια ομάδα αρχείων. Από αυτήν μερικά αρχεία είναι απαραίτητα και κάποια άλλα είναι προεραιτικά.

Τα απαραίτητα αρχεία για τον ορισμό του shapefile είναι τα .shp, .shx και .dbf. Τα δυο πρώτα σχετίζονται με την γεωγραφική πληροφορία του αρχείου ενώ το τελευταίο με την πληροφορία την οποία απεικονίζει. Για παράδειγμα, ένα σημειακό shapefile που απεικονίζει σημεία ενδιαφέροντος αποθηκεύει το είδος του αρχείου (σημειακό) και τις συν/νες των σημείων στα .shp και .shx. Αν τώρα για κάθε σημείο υπάρχει ένας μοναδικός κωδικός, ένας κωδικός για τον δήμο που ανήκει και μια περιγραφή σχετικά με το τι σημείο είναι (οργανισμός, μουσείο, κινηματογράφος κ.λπ.) αυτά αποθηκεύονται στο .dbf, το οποίο ουσιαστικά είναι ένα αρχείο-πίνακας με τόσες σειρές όσα τα σημεία και τόσες στήλες όσα τα περιγραφικά χαρακτηριστικά τους (σε αυτή την περίπτωση 3).

Για την περίπτωση των γραμμικών και επιφανειακών στοιχείων στο .shp και .shx αποθηκεύονται όλες οι συν/νες της polyline που ορίζει την κάθε γραμμή.

Από τα προεραιτικά αρχεία ιδιαίτερο ενδιαφέρον παρουσιάζει το αρχείο .prj (από το projection), το οποίο περιέχει την πληροφορία σχετικά με το σύστημα συντεταγμένων του αρχείου.

Το ogrinfo

Πριν προχωρήσουμε παρακάτω καλό είναι να αναφερθεί το πρόγραμμα ogrinfo, το οποίο δίνει βασικές πληροφορίες για ένα shapefile. Από ένα τερματικό (Applications → Accessories → Terminal) με την εντολή:

  ogrinfo -so /path/to/the/shapefile.shp shapefile

Για παράδειγμα το αρχείο roads.shp με την εντολή

 ogrinfo -so /home/user/roads.shp roads

επιστρέφει:

  INFO: Open of `/home/vagvaf/Data/greece.shp/roads.shp'
     using driver `ESRI Shapefile' successful.
  Layer name: roads
  Geometry: Line String
  Feature Count: 88925
  Extent: (19.632112, 34.929228) - (27.353774, 41.754906)
  Layer SRS WKT:
  GEOGCS["GCS_WGS_1984",
       DATUM["WGS_1984",
            SPHEROID["WGS_1984",6378137,298.257223563]],
       PRIMEM["Greenwich",0],
       UNIT["Degree",0.017453292519943295]]
  osm_id: Real (11.0)
  name: String (48.0)
  ref: String (16.0)
  type: String (16.0)
  oneway: Integer (1.0)
  maxspeed: Integer (3.0)
 
geo-processing.txt · Τελευταία τροποποίηση: 2014/05/07 14:25 από Sairin_Lote