Αυτοματοποιημένη χαρτογράφηση της σεισμικής δραστηριότητας με βιβλιοθήκες ΕΛΛΑΚ.

Αρκετοί χρήστες, κυρίως αρχάριοι, αναρωτιούνται γιατί να χρησιμοποιήσει κάποιος μια γλώσσα προγραμματισμού και να συντάξει ένα σενάριο (script), πολύπλοκων αρκετές φορές, εντολών για να παράγει ένα χάρτη από την στιγμή που υπάρχει λογισμικό με γραφικό περιβάλλον που με ευκολία μπορεί να οδηγήσει στο ίδιο αποτέλεσμα; Η απάντηση είναι, μεταξύ άλλων, γιατί θέλει να αυτοματοποιήσει την διαδικασία και να αποφύγει ο χρήστης την επανάληψη εργασιών που πολλές φορές στοιχίζει σε χρόνο και χρήμα. Με την σύνταξη του κατάλληλου σεναρίου ο υπολογιστής αναλαμβάνει να εκτελέσει όλη την εργασία και να την επαναλαμβάνει όποτε ο χρήστης του ορίσει.

Για να γίνει πιο κατανοητή η δύναμη της αυτοματοποιημένης χαρτογραφίας παρατίθεται ένα παράδειγμα χαρτογράφησης των σεισμών σε παγκόσμιο επίπεδο.

Τα δεδομένα για το υπόβαθρο του χάρτη προέρχονται από την ιστοσελίδα Natural Earth ενώ τα δεδομένα με τα επίκεντρα των σεισμών προέρχονται από το United States Geological Survey (USGS).

Διαλέξαμε το αρχείο που περιλαμβάνει τους σεισμούς των τελευταίων 7 ημερών μεγέθους άνω των 2,5 Richter.

Σκοπός του σεναρίου (script) είναι:

• Να κατεβάσει από το διαδίκτυο τα απαραίτητα δεδομένα σε μορφή CSV.
• Να μετατρέψει το αρχείο CSV σε Shapefile με την βιβλιοθήκη ogr.
• Να παράγει τον τελικό χάρτη που θα συνδυάζει το γεωφυσικό υπόβαθρο με τα επικεντρα των σεισμών με βοήθεια της βιβλιοθήκης Mapnik.
• Να φορτώσει τον χάρτη σε μορφή PNG στην τρέχουσα ιστοσελίδα.
• Να εκτελείται επαναλαμβανόμενα σε καθημερινή βάση ώστε ο χάρτης να ανεβαίνει ενημερωμένος στην σελίδα.
• Το σενάριο συντάχθηκε με την Python στο λειτουργικό σύστημα Ubuntu Linux 11.04 (Natty Narwhal).

Το σενάριο:

# -*- coding: utf-8 -*-
#!/usr/bin/python
 
#Εισάγουμε τις απαραίτητες βιβλιοθήκες
import os,  urllib, mapnik,  Image,  csv,  ftplib,  datetime
from osgeo import ogr
from osgeo import osr
from StringIO import StringIO
 
import PIL
from PIL import ImageFont
from PIL import Image
from PIL import ImageDraw
 
#Κατεβάζουμε τα απαραίτητα δεδομένα σε csv μορφή
print "Downloading data file..."
now = datetime.datetime.now()
url ='http://earthquake.usgs.gov/earthquakes/catalogs/eqs7day-M2.5.txt'
file = 'earthquakes.txt'
urllib.urlretrieve( url,file)
 
#Δημιουργία του shapefile
outShapeFile = 'earthquakes.shp'
drv = ogr.GetDriverByName('ESRI Shapefile')
 
if os.path.exists(str(outShapeFile)):
drv.DeleteDataSource(str(outShapeFile))
print "Deleting shapefile that already exists: " + outShapeFile
 
# Δημιουργία SpatialReference
t_srs = osr.SpatialReference()
t_srs.SetFromUserInput('WGS84')
 
ds = drv.CreateDataSource(outShapeFile)
print "Creating shapefile: " + outShapeFile
layer = ds.CreateLayer(ds.GetName(), geom_type = ogr.wkbPoint, srs = t_srs)
 
#Τα πεδία με το όνομα τους και τις ιδιότητές τους
fields=[['Src',  ogr.OFTString, 5],
['Eqid',  ogr.OFTString , 15],
['Version',  ogr.OFTInteger] ,
['Datetime',  ogr.OFTString, 255] ,
['Lat',  ogr.OFTReal],
['Lon',  ogr.OFTReal],
['Magnitude',  ogr.OFTReal],
['Depth',  ogr.OFTReal],
['NST',  ogr.OFTReal],
['Region',  ogr.OFTString, 100]
]
 
#Δημιουργία πεδίων
for field in fields:
fieldname=field[0]
fieldproperty=field[1]
print "Creating Field "  + fieldname
ofield  = ogr.FieldDefn(fieldname, fieldproperty)
if fieldproperty ==4:#Ειδικό property αν είναι string
fieldlen=field[2]
ofield.SetWidth(fieldlen)
layer.CreateField(ofield)
 
print "Creating geometry object..."
geom = ogr.Geometry(type=ogr.wkbPoint)
 
# αρχικοποίηση μιας απαραίτητης μεταβλητής
Line = ""
 
#Ανάγνωση αρχείου csv
f=open(str('earthquakes.txt'), 'r')
f.readline()
while Line.isalnum:
Line = f.readline()
if not Line: break
data = StringIO(Line)
reader = csv.reader(data, delimiter=',')
for row in reader:
geom.AddPoint(float(row[5]),float(row[4]))#Προσθήκη ενός σημείου
feat = ogr.Feature(feature_def=layer.GetLayerDefn())
print "Importing shape into feature."
feat.SetGeometry(geom)
#Προσθήκη δεδομένων στο feature
print "Importing attributes into feature."
feat.SetField('Src',row[0] )
feat.SetField('Eqid',row[1] )
feat.SetField('Version',row[2] )
feat.SetField('Datetime',row[3] )
feat.SetField('Lat',row[4] )
feat.SetField('Lon',row[5] )
feat.SetField('Magnitude',row[6] )
feat.SetField('Depth',row[7] )
feat.SetField('NST',row[8] )
feat.SetField('Region',row[9] )
layer.CreateFeature(feat)
 
#Διαγραφή objects
feat.Destroy()
ds.Destroy()
 
#Δημιουργία του χάρτη με το mapnik
 
mapfile = 'world_styles.xml' #απαραίτητο xml με τις πηγές των αρχείων και διάφορες ιδιότητές τους
map_output = 'earthquakes.png'
m = mapnik.Map(1024, 550)
mapnik.load_map(m, mapfile)
bbox = mapnik.Envelope(mapnik.Coord(-180.0, -90.0), mapnik.Coord(180.0, 90.0))
m.zoom_to_box(bbox)
mapnik.render_to_file(m, map_output,  "png") #εξαγωγή σε png
 
#Αναγραφή της ημερομηνίας ενημέρωσης του χάρτη.
font = ImageFont.truetype("/usr/share/fonts/truetype/ttf-dejavu/DejaVuSansCondensed.ttf",12)
infile = "earthquakes.png"
img=Image.open(infile)
draw = ImageDraw.Draw(img)
 
curdate= now.strftime("%d-%m-%Y %H:%M")
draw.text((10, 0),u"Τελευταία ενημέρωση:" + curdate + ", by geographer.gr",(0, 0,0),font=font)
draw = ImageDraw.Draw(img)
img.save("earthquakes.png")
 
#Δημιουργία thumbnail
print "Generate the thumbnail"
thumbfile="earthquakes_404x196_a0da98766e696ec356a29088bc2d8e16.jpg"
image = Image.open(map_output)
image = image.resize((500, 250), Image.ANTIALIAS)
image.save(thumbfile)
 
#Ανέβασμα των χαρτών με ftp στους απαραίτητους φακέλους στο geographer.gr
print "Upload data on geographer.gr"
ftp = ftplib.FTP("geographer.gr")
print ftp.login("myusername", "mypass")
print ftp.cwd('/public_html/images/stories/earthquakes')
print ftp.storbinary("STOR " + map_output, open(map_output, "rb"))
 
print ftp.cwd('/public_html/plugins/content/fboxbot/thumbs')
print ftp.storbinary("STOR " + thumbfile, open(thumbfile, "rb"))
print 'Closing FTP connection'
print ftp.close()
 
print "Job done!"

Όπως φαίνεται παραπάνω το mapnik χρησιμοποιεί για την δημιουργία του χάρτη ένα αρχείο xml όπου περιέχει την προέλευση των αρχείων των δεδομένων και διάφορες ιδιότητες τους όπως το γεωγραφικό σύστημα αναφορά, σύμβολα και κλάσεις. Το περιεχόμενο του αρχείου δίδεται παρακάτω:

<?xml version="1.0" encoding="utf-8"?>
<Map srs="+init=epsg:4326" bgcolor="rgb(255,255,255)">
<Style name="NE2_50M_SR_W_style">
<Rule>
<RasterSymbolizer>
<CssParameter name="scaling">bilinear</CssParameter>
</RasterSymbolizer>
</Rule>
</Style>
 
<Style name="earthquakes_style">
<Rule>
<Filter>(([Magnitude]&amp;gt;=2.5) and ([Magnitude]&amp;lt;=3.5))</Filter>
<PointSymbolizer allow_overlap="yes" file="symbols/sym_0.png" type="png" width="7" height="7" opacity="1"></PointSymbolizer>
</Rule>
<Rule>
<Filter>(([Magnitude]&amp;gt;3.5) and ([Magnitude]&amp;lt;=4.5))</Filter>
<PointSymbolizer allow_overlap="yes" file="symbols/sym_1.png" type="png" width="9" height="9" opacity="1"></PointSymbolizer>
</Rule>
<Rule>
<Filter>(([Magnitude]&amp;gt;4.5) and ([Magnitude]&amp;lt;=5))</Filter>
<PointSymbolizer allow_overlap="yes" file="symbols/sym_2.png" type="png" width="13" height="13" opacity="1"></PointSymbolizer>
</Rule>
<Rule>
<Filter>(([Magnitude]&amp;gt;5) and ([Magnitude]&amp;lt;=5.5))</Filter>
<PointSymbolizer allow_overlap="yes" file="symbols/sym_3.png" type="png" width="15" height="15" opacity="1"></PointSymbolizer>
</Rule>
<Rule>
<Filter>(([Magnitude]&amp;gt;5.5) and ([Magnitude]&amp;lt;=6.2))</Filter>
<PointSymbolizer allow_overlap="yes" file="symbols/sym_4.png" type="png" width="21" height="21" opacity="1"></PointSymbolizer>
</Rule>
</Style>
 
<Layer name="HYP_50M_SR_W" srs="+init=epsg:4326">
<StyleName>NE2_50M_SR_W_style</StyleName>
<Datasource>
<Parameter name="file">NE2_50M_SR_W/NE2_50M_SR_W.tif</Parameter>
<Parameter name="type">gdal</Parameter>
</Datasource>
</Layer>
<Layer name="earthquakes" srs="+init=epsg:4326">
<StyleName>earthquakes_style</StyleName>
<Datasource>
<Parameter name="file">earthquakes.shp</Parameter>
<Parameter name="type">shape</Parameter>
</Datasource>
</Layer>
 
</Map>
 
<!-- nik2img.py mapnik.xml out.png -d 1017 308 -e -302.42451763 -76.2454581721 301.57881763 105.848158172 -->

Για την καθημερινή επαναλαμβανόμενη εκτέλεση του σεναρίου (π.χ. στις 8:50 κάθε πρωί) χρησιμοποιούμε την υπηρεσία cron του linux, όπου εισάγουμε την παρακάτω εγγραφή στο αρχείο crontab:

50 8 * * * cd ~/GIS/earthquakes;/usr/bin/python