Jython et le géospatial: exemples avec GeoScript (GeoTools) , JTS ( Java Topology Suite), OpenJump, gvSIG, WorldWind ou Geoserver

lun 03-12-2012 Martin Laloux

SIG OpenSource

python

Tout le monde connait le Python standard (ou « C Python »), qui s'exécute dans un environnement C, mais, peu connaissent ou utilisent ses autres implémentations (voir « Python: géospatial, dialectes (standard, pour ESRI, pour FME, pour GvSIG etc... ») .

Parmi celles-ci, il y a Jython, qui est son implémentation en Java et qui s'exécute dans une Machine Virtuelle Java (JVM). Dans le monde géospatial, il est utilisé comme langage de script pour les SIGs écrits en Java, comme gvSIG ou OpenJump et pour un serveur cartographique, GeoServer.

Mais, tout comme Python, Jython peut-être utilisé directement pour traiter les données géospatiales. Il peut utiliser tous les modules Python écrits en Python pur, mais, comme déjà signalé, il sera impossible d'utiliser des modules écrits en C (GDAL, numpy,...). En contrepartie, Jython pourra utiliser toutes les classes des librairies Java existantes. Ses caractéristiques sont les suivantes (fr.wikipedia.org/wiki/Jython):

le code Python est compilé en bytecode Java ;
les classes Python peuvent hériter des classes Java;
possibilité d'exécuter des codes Python durant le fonctionnement d’un programme Java (utilisée par GvSIG ou OpenJump) ;
utilisation d’objets Java dans le code Python.

Dans Jython, la variable CLASSPATH est fondamentale, car c'est elle qui permet l'utilisation de n'importe quelle classe Java. Ainsi si je veux utiliser le driver JDBC natif (postgresql-xxx.jdbc3.jar) pour me connecter à PostgreSQL (voir www.chicoree.fr/w/Utiliser_JDBC_%C3%A0_partir_de_Jython):

connexion

# ajustement du CLASSPATH de Jython pour l'importation des classes du driver .jar
$ export CLASSPATH="/Library/Java/Extensions/postgresql-9.0-801.jdbc3.jar"
 
$ jython
# importation des classes (première utilisation, après cela ne sera plus signalé)
*sys-package-mgr*: processing new jar, '/Library/Java/Extensions/postgresql-9.0-801.jdbc3.jar'
 
Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06) 
[Java HotSpot(TM) 64-Bit Server VM (Apple Inc.)] on java1.6.0_24
Type "help", "copyright", "credits" or "license" for more information.
>>> from org.postgresql import Driver
>>> 
# si aucune erreur n'est signalée, c'est que ça marche

3303 lectures

Je n'irai pas plus loin dans la présentation de Jython (installation sur les diverses plates-formes , principes, etc.) renvoyant les lecteurs à des sites plus spécialisés comme opikanoba.org/java/java-et-python (en français, mais, il y en a beaucoup d'autres). Pour débuter, il y a un livre présenté en Open Source avec une licence Creative Commons 3.0. , « The Definitive Guide to Jython » (www.jython.org/jythonbook/en/1.0/), tous les livres et articles présentés à wiki.python.org/jython/JythonBibliography, des blogs ou des forums comme le Forum des développeurs (www.developpez.net/forums/).

Je me limiterai ici à son application dans le domaine géospatial et dans la suite, nous traiterons:

l'utilisation de Jython seul avec des modules géospatiaux (GeoScript, JTS);

l'utilisation de Jython comme langage de script de logiciels SIG (OpenJump, GvSIG, Worldwind;

l'utilisation de Jython avec un serveur cartographique (GeoServer)

Utilisation de Jython seul

GeoScript ou l'équivalent de Shapely pour Jython

Shapely, dépendant de la librairie C GEOS, ne peut donc pas être utilisé, mais, il a son équivalent pour Jython, GeoScript, proposé par GeoTools « The Open Source Java GIS Toolkit » qui est un produit de l'OSGeo. GeoScript permet d'utiliser les classes de GeoTools en Jython, JavaScript, Scala et Groovy.

Comment l'installer:

La manière la plus simple est d'installer easy_install (c'est possible depuis les versions 2.5.x de Jython). Pour ce faire, il faut télécharger le script peak.telecommunity.com/dist/ez_setup.py et ensuite:

jython ez_setup.py

Ensuite, après téléchargement et décompression de GeoScript (geoscript.org/py/quickstart.html#quickstart)

jython setup.py install (easy_install va vous installer toutes les dépendances)

Ce qui va vous créer (entre autres choses) deux fichier dans le dossier bin du dossier Jython2.5.x, geoscript et geoscript-classpath

Utilisation

Il y a 2 manières d'utiliser GeoScript:

soit à partir de l'application geoscript

avec geoscript

$ geoscript
Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06) 
[Java HotSpot(TM) 64-Bit Server VM (Apple Inc.)] on java1.6.0_24
Type "help", "copyright", "credits" or "license" for more information.
>>>

4022 lectures

soit à partir de jython lui-même

avec jython

# ajustement du CLASSPATH
export CLASSPATH=`/.../jython2.5.2/bin/geoscript-classpath`
$ jython
Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06) 
[Java HotSpot(TM) 64-Bit Server VM (Apple Inc.)] on java1.6.0_24
Type "help", "copyright", "credits" or "license" for more information.
>>> import geoscript
>>> 
 
# si vous ne le faites pas
$ jython
Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06) 
[Java HotSpot(TM) 64-Bit Server VM (Apple Inc.)] on java1.6.0_24
Type "help", "copyright", "credits" or "license" for more information.
>>> import geoscript
Error: Could not find GeoTools libraries on classpath.

4261 lectures

Ensuite, c'est pratiquement comme avec le module Shapely (voir www.portailsig.org/content/python-le-module-shapely-geometries-predicats-spatiaux-analyse-spatiale-matrices-de-clementi)

traitement

>>> from geoscript import geom
>>> point = geom.readWKT('POINT(45 20)')
>>> point
POINT (45 20)
# ou
>>> point = geom.Point(4.607817,50.664930)
>>> point
POINT (4.607817 50.66493)
>>> # ou
>>> wkt = str(point)
>>> wkt
'POINT (4.607817 50.66493)'

3455 lectures

Mais, en plus, GeoScript permet:

les changements de système de projection:

projections

>>> from geoscript import proj
>>> point = geom.Point(4.607817,50.664930)               
>>> point2 = proj.transform(point, 'epsg:4326', 'epsg:31370')
>>> point2
POINT (166901.5939415752 150507.1768997945)

3757 lectures

l'importation directe des shapefiles et leur création (et des couches PostGIS):

Shapefiles

>>> from geoscript.layer import Shapefile
>>> macouche = Shapefile('.../montest.shp')
>>> from geoscript.layer import Shapefile
>>> pntGeoms = [f.geom for f in macouche.features()]
>>> pntGeoms
[MULTILINESTRING ((0.0316820276497696 0.6129032258064516, 0.1157834101382488 0.7177419354838709, 0.2217741935483871 0.683179723502304, 0.2839861751152074 0.5656682027649769, 0.3680875576036867 0.6048387096774194, 0.4326036866359447 0.7085253456221199, 0.5293778801843319 0.6071428571428571)), MULTILINESTRING ((0.3162442396313365 0.8940092165898617, 0.2321428571428572 0.8421658986175116, 0.2459677419354839 0.7039170506912441, 0.1756912442396313 0.5794930875576036, 0.2736175115207373 0.4331797235023041))]

4169 lectures

leur représentation graphique (directement par Java):
plot
>>> from geoscript.viewer import plot
>>> plot(pntGeoms)
3777 lectures
etc. (voir geoscript.org/py/index.html ou geoscriptblog.blogspot.com/ pour de nombreux exemples)

avec la librairie Java, JTS - Java Topology Suite

La Java Topology Suite est la librairie Java pour le traitement des géométries en 2 dimensions. La librairie C GEOS en est dérivée. Elle est utilisée par GeoTools, GeoServer, tous les Jump's, GvSIG ou Udig par exemple (voir tsusiatsoftware.net/jts/jts-links.html#systems). GeoScript est basée sur cette librairie, mais toutes ses classes n'ont pas encore été portées.

JTS ne possède pas d'implémentation Python, mais nous avons vu que Jython permet d'employer n'importe quelle classe Java. Alors certains ont agi...

Jython et JTS

>>> import sys
>>> # autre manière d'importer des classes (ici JTS)
>>> sys.path.append('.../jts-1.8.0/lib/jts-1.8.jar')
>>> # test d'une classe de JTS
>>> from com.vividsolutions.jts.geom import Point
>>>

4526 lectures

Le problème est de trouver comment appeler les classes Java de JTS en Jython. Cela est réalisé par certains à coup d'essais et d'erreurs. Il est cependant nécessaire de bien connaitre la bibliothèque Java au préalable. C'est le cas de Carl Trachte sur son site geoscriptblog.blogspot.com/2010/06/unwrapped-jts-with-python.html, qui offre de très belles applications.

résultat d'une triangulation de Delaunay avec JTS

Utilisation de Jython comme langage de script de logiciels SIG

Jython peut aussi, en théorie, être utilisé comme langage de script dans des applications Java comme gvSig, OpenJump (et la grande famille des Jump's, voir www.portailsig.org/content/sortie-de-la-version-151-de-openjump-gis-et-catalogue-de-tous-les-logiciels-derives-de-jump) ou Worldwind Java. Pourquoi en théorie ? Parce que son utilisation nécessite une bonne connaissance des classes Java de l'application et donc de Java et qu'il y a très peu d'exemples disponibles.

avec OpenJump (et la famille des Jump's)

Nous commençons par OpenJump car son module Jython est l'application directe de l'utilisation de la Java Topology Suite (JTS).

Malheureusement, bien que très puissant, il faut « chipoter » pour trouver comment ce module fonctionne. Il n'y a ni tutoriel ni exemples. J'y suis arrivé en analysant systématiquement les fonctions des objets obtenus et la connaissance de la manière de fonctionner des modules Python dans le domaine geospatial. Comme cet aspect est mal connu pour OpenJump, une petite synthèse des résultats obtenus:

traitements avec le module Jython d'OpenJump

# importation du module principal d'OpenJump
from org.openjump.util.python.JUMP_GIS_Framework import *
# obtention des couches sélectionnées
couches = getSelectedLayers()
# fonctions disponibles pour les couches sélectionnées:
dir(couches)
['add', 'addAll', 'class', 'clear', 'clone', 'contains', 'containsAll', 'empty', 'ensureCapacity', 'equals', 'get', 'getClass', 'hashCode', 'indexOf', 'isEmpty', 'iterator', 'lastIndexOf', 'listIterator', 'notify', 'notifyAll', 'remove', 'removeAll', 'retainAll', 'set', 'size', 'subList', 'toArray', 'toString', 'trimToSize', 'wait']
 
# itération à travers les couches sélectionnées
for i in couches.iterator():
    print couche.name
    ...
# fonctions disponibles pour une couche:	
dir(couche)
['FIRING_APPEARANCE_CHANGED_ON_ATTRIBUTE_CHANGE', 'addStyle', 'addUndo', 'basicStyle', 'blackboard', 'class', 'cloneStyles', 'dataSourceQuery', 'defaultLineColor', 'description', 'dispose', 'drawingLast', 'editable', 'equals', 'featureCollection', 'featureCollectionModified', 'featureCollectionWrapper', 'fireAppearanceChanged', 'fireLayerChanged', 'getBasicStyle', 'getBlackboard', 'getClass', 'getDataSourceQuery', 'getDescription', 'getFeatureCollectionWrapper', 'getLabelStyle', 'getLayerManager', 'getMaxScale', 'getMinScale', 'getName', 'getStyle', 'getStyles', 'getTask', 'getVertexStyle', 'hasReadableDataSource', 'hashCode', 'isDrawingLast', 'isEditable', 'isFeatureCollectionModified', 'isReadonly', 'isScaleDependentRenderingEnabled', 'isSelectable', 'isSynchronizingLineColor', 'isVisible', 'labelStyle', 'layerManager', 'maxScale', 'minScale', 'name', 'notify', 'notifyAll', 'readonly', 'removeStyle', 'scaleDependentRenderingEnabled', 'selectable', 'setDataSourceQuery', 'setDescription', 'setDrawingLast', 'setEditable', 'setFeatureCollection', 'setFeatureCollectionModified', 'setLayerManager', 'setMaxScale', 'setMinScale', 'setName', 'setReadonly', 'setScaleDependentRenderingEnabled', 'setSelectable', 'setStyles', 'setSynchronizingLineColor', 'setVisible', 'styles', 'synchronizingLineColor', 'task', 'toString', 'tryToInvalidateEnvelope', 'vertexStyle', 'visible', 'wait']
 
### Applications ###
# itération à travers les géométries d'une couche
fc = couches[0].featureCollectionWrapper
for elem in range(fc.size):
    obj = fc.features[elem]	
	geom = obj.getGeometry()
	...	
 
# fonctions disponibles pour la géométrie d'un élément (la 3D est reconnue par défaut):
dir(geom)
['SRID', 'apply', 'area', 'boundary', 'boundaryDimension', 'buffer', 'centroid', 'class', 'clone', 'compareTo', 'contains', 'convexHull', 'coordinate', 'coordinates', 'coveredBy', 'covers', 'crosses', 'difference', 'dimension', 'disjoint', 'distance', 'empty', 'envelope', 'envelopeInternal', 'equals', 'equalsExact', 'equalsNorm', 'equalsTopo', 'exteriorRing', 'factory', 'geometryChanged', 'geometryType', 'getArea', 'getBoundary', 'getBoundaryDimension', 'getCentroid', 'getClass', 'getCoordinate', 'getCoordinates', 'getDimension', 'getEnvelope', 'getEnvelopeInternal', 'getExteriorRing', 'getFactory', 'getGeometryN', 'getGeometryType', 'getInteriorPoint', 'getInteriorRingN', 'getLength', 'getNumGeometries', 'getNumInteriorRing', 'getNumPoints', 'getPrecisionModel', 'getSRID', 'getUserData', 'hashCode', 'interiorPoint', 'intersection', 'intersects', 'isEmpty', 'isRectangle', 'isSimple', 'isValid', 'isWithinDistance', 'length', 'norm', 'normalize', 'notify', 'notifyAll', 'numGeometries', 'numInteriorRing', 'numPoints', 'overlaps', 'precisionModel', 'rectangle', 'relate', 'reverse', 'setSRID', 'setUserData', 'simple', 'symDifference', 'toString', 'toText', 'touches', 'union', 'userData', 'valid', 'wait', 'within']
# d'où :
geom.geometryType
'Polygon'
# la géométrie
geom.envelope
POLYGON ((226023.7 18173, 226023.7 28173.5, 234026.3 28173.5, 234026.3 18173, 226023.7 18173))
# nombre de points
geom.getNumPoints()
5
# les coordonnées des points de la géométrie 
geom.coordinates
array(com.vividsolutions.jts.geom.Coordinate,[(226023.7, 28173.5, NaN), (234025.0, 28173.5, NaN), (234026.3, 18173.1, NaN), (226024.5, 18173.0, NaN), (226023.7, 28173.5, NaN)]) 
geom.coordinate
(226023.7, 28173.5, NaN)
# le centroide
geom.getCentroid()
POINT (230024.86833801787 23173.222923277397)
# la limite extérieure de la géométrie
geom.getBoundaryDimension()
1
geom.getBoundary()
LINEARRING (226023.7 28173.5, 234025 28173.5, 234026.3 18173.1, 226024.5 18173, 226023.7 28173.5)
# ou
geom.exteriorRing
LINEARRING (226023.7 28173.5, 234025 28173.5, 234026.3 18173.1, 226024.5 18173, 226023.7 28173.5)
 
# les attributs
fonctions disponibles pour les attributs d'un élément:
>>> dir(obj)
['ID', 'attributes', 'class', 'clone', 'compare', 'compareTo', 'equals', 'geometry', 'getAttribute', 'getAttributes', 'getClass', 'getDouble', 'getGeometry', 'getID', 'getInteger', 'getSchema', 'getString', 'hashCode', 'isModified', 'modified', 'notify', 'notifyAll', 'schema', 'setAttribute', 'setAttributes', 'setGeometry', 'setModified', 'setSchema', 'toString', 'wait']
 
# d'où
obj.attributes
array(java.lang.Object,[POLYGON ((226023.7 28173.5, 234025 28173.5, 234026.3 18173.1, 226024.5 18173, 226023.7 28173.5)), '71/5', '71', '225W', 'LAMORTEAU', 'LAMORTEAU - RUETTE', '71/5-6']) 
# on retrouve donc la géométrie
obj.getGeometry()
POLYGON ((226023.7 28173.5, 234025 28173.5, 234026.3 18173.1, 226024.5 18173, 226023.7 28173.5))
# ou
obj.getAttribute(0)
POLYGON ((226023.7 28173.5, 234025 28173.5, 234026.3 18173.1, 226024.5 18173, 226023.7 28173.5))
# attributs éventuels (fonction de la liste obtenue)
obj.getAttribute(1)
'71/5'
obj.getString(3)
'225W'

2997 lectures

Et donc, il y a moyen d'utiliser directement toutes les fonctions de la JTS:

Il est aussi possible de remarquer que bien qu'OpenJump ne gère pas les projections, son module Jython permet théoriquement de le faire (SRID).

Il a aussi moyen de créer des extensions/modules en Jython (voir mekandizim.mersin.edu.tr/belgeler/SSA-OJ.pdf, par exemple, et son module disponible à mekandizim.mersin.edu.tr/eklentiler/eklentiler.html)

avec gvSIG:

Le traitement est comparable à celui d'OpenJump mais il sera moins développé ici du fait d'un changement de module avec la future version 2.0 en développement.

avec les versions 1.x, c'est le module gvsiglib. Il faut néanmoins fouiller dans les listes de gvSIG de l'OSGeo (osgeo-org.1803224.n2.nabble.com/gvSIG-project-f2000051.html) pour trouver des exemples. Il y a un « scripting tutorial » ( forge.osor.eu/docman/view.php/89/359/gvsig-1_0-scripting-tutorial-v2-en.zip), très léger, ou des blogs comme dominoc925.blogspot.be/2011/07/simple-gvsig-jython-console-script-to.html.

avec la future version 2.0, c'est le nouveau module gvsig, beaucoup plus « pythonesque » (traduction libre de "pythonic") mais il faut, pour l'utiliser, télécharger la dernière version disponible (www.gvsig.org/web/projects/gvsig-desktop/devel/gvsig/2_0_0/version) et installer l'extension Scripting Framework. Il y a des versions pour Linux et Windows, mais rien pour Mac OS X pour l'instant. Patience, Agustín Diez Castillo (gvsigmac.blogspot.be/ ) y travaille (communication personnelle).
Exemple d'utilisation reprise du blog de gvSIG (blog.gvsig.org/2012/10/22/scripting-exprime-tu-gvsig/):

from gvsig import *
from commonsdialog import *
def main():
  layer = currentLayer()
  if layer == None:
    msgbox("Vous devez choisir une couche!", "AVIS", 1)
    return
  emax = 0.0
  emin = 0.0
  for feature in layer.features():
    if feature.ELEVATION &gt; emax :
      emax = feature.ELEVATION
    if feature.ELEVATION &lt; emin or emin ==0.0:
      emin = feature.ELEVATION
  msgbox("Altitude maximum=%s et minimum=%s" % (emax, emin), "Elevation", 0)

Résultat:

image reprise de blog.gvsig.org/2012/10/22/scripting-exprime-tu-gvsig/

Il est aussi possible de créer des extensions/modules en Jython (voir mekandizim.mersin.edu.tr/belgeler/book.pdf et la même référence que pour OpenJump)

avec Worldwind Java

Je n'ai trouvé qu'un seul exemple d'utilisation à www.perrygeo.net/wordpress/.

Utlisation de Jython avec un serveur cartographique: GeoServer

GeoServer a aussi lancé son module Jython (geoserver.org/display/GEOS/Python+Scripting+Extension). Encore expérimental, il permet d'interroger Geoserver et même d'interagir avec lui.

from geoserver import Catalog
cat = Catalog('aff_geol')
print cat.stores()
[u'affleurements', u'cartes']
st = cat['affleurements']
print st.layers()
[u'affs']
l = st['affs']
print l.count()
3214
print l.bounds()
(107206.74,148102.45,114110.46,157514.40, EPSG:31370)

Conclusions

Tant qu'il s'agit de modules « clés en main », comme GeoScript, Jython est facile à utiliser, mais dès que l'on veut utiliser directement une librairie Java, une bonne connaissance de ce langage s'impose. Le même constat s'adresse aussi aux langages de scripts des SIGs.

Le nouveau module gvsig de la version 2 de gvSIG parait néanmoins se rapprocher beaucoup plus d'une version Python (sans devoir connaître le détail des classes Java dont il dépend). Et finalement j'ai été fort impressionné par les possibilités du module d'OpenJump (logiciel à mon avis sous-évalué).

Je regrette toujours qu'il ne soit pas possible d'utiliser tout en même temps, à savoir PyQGIS dans gvSIG et l'inverse ou les deux dans un shell, par exemple....

Mais, en pratique, connaître Python + Java, le tout sans un tutoriel détaillé, rend son utilisation quelque peu difficile pour le commun des mortels ...

Tous les traitements ont été effectués sur Mac OS X avec Jython 2.5.2, GeoScript version 1.1, OpenJump version 1.5.1 rev.2721 et diverses implémentations de gvSIG.