Introduction

Dans mon précédent article sur QGis et Oracle Spatial, j'avais indiqué comment créer une couche dans Oracle pour la rendre "compatible" avec QGis. Si vous avez de nombreuses couches à importer, écrire manuellement ce code SQL peut entraîner des erreurs. De plus, certains paramètres tels que l'emprise de la bounding-box sont fastidieux à saisir et à calculer.

J'ai donc imaginé un moyen de rendre cet import plus facile. Ce travail est sans aucun doute facilement améliorable mais il a l'intérêt de proposer une méthode complètement intégrée à QGis sans devoir développer un outil tiers ou un plugin dédié.

La technique

Depuis la version 2.0 de QGis, il existe un plugin un peu particulier. C'est le plugin Processing, auparavant appelé Sextante. Il propose de reproduire le concept des géotraitements qu'on trouve dans ESRI ArcGIS Desktop. La ressemblance est d'ailleurs assez frappante lorsqu'on ouvre les deux interfaces graphiques pour comparer.

Processing met à disposition un ensemble de scripts Python qui forment des GéoAlgorithmes. Ces derniers réalisent des traitements sur des couches en entrée et envoient de la donnée en sortie, dans une autre couche ou dans d'autres types de sortie (fichier à plat, variable texte ou nombre, etc.). Il est possible de coupler ces GéoAlgorithmes entre eux pour former des files de traitement. Processing propose d'ailleurs une interface graphique dédiée qui vous permet de composer une ou plusieurs chaînes de traitement.

Composeur de chaîne de traitement

Processing propose également des scripts en Python qui sont exécutés avec une boîte de dialogue simple. Pour configurer cette boîte de dialogue, il suffit d'ajouter des commentaires spéciaux dans le code (marqués avec un double #).

Installer le script

Pour installer le script, le plus simple est de copier le fichier dans le répertoire utilisateur dédié à Processing:

  • ~/.qgis2/processing/scripts/ sous GNU/Linux
  • C:\Users\UTILISATEUR.qgis2\processing\scripts\ sous MS-Windows

A la suite de cet ajout et d'un redémarrage de QGis, le script devrait être disponible dans la fenêtre des géotraitements:

Fenêtre Processing

Le script

Voici le code du script en Python 2.7, version actuellement utilisée par la version 2.2 de QGis. Vous pourrez noter que les paramètres de la boîte de dialogue sont décrits dès le début du script.

# -*- coding: utf-8 -*-

##[Oracle]=group
##input=vector
##table_name=string
##schema=string
##output=output file D:/TEMP/test.sql
##drop_table=boolean False
##add_spatial_index=boolean False
##primary_key=field input
##uppercase_field_names=boolean False
##srid=number 2154
##geometry_type=selection UNKNOWN;POINT;MULTIPOINT;POINT+MULTIPOINT;LINESTRING;MULTILINESTRING;LINESTRING+MULTILINESTRING;POLYGON;MULTIPOLYGON;POLYGON+MULTIPOLYGON;Scan the layer to determine
##tolerance=number 0.01
##inject_data=boolean False


# This script generates an SQL File to create a table in Oracle Spatial
# You can also import the data

from qgis.core import *
from qgis.core import QgsCoordinateTransform
from PyQt4.QtCore import *
from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException

# Get some infos from the layer
layer = processing.getObject(input)
provider = layer.dataProvider()
fields = provider.fields()

# Get the extent and manage the SRID transformations:
crsSrc = layer.crs()
crsDest = QgsCoordinateReferenceSystem(int(srid), QgsCoordinateReferenceSystem.EpsgCrsId)
crsTransform = QgsCoordinateTransform(crsSrc, crsDest)

lower_bound = crsTransform.transform(layer.extent().xMinimum(), layer.extent().yMinimum())
upper_bound = crsTransform.transform(layer.extent().xMaximum(),layer.extent().yMaximum())

extent = { 'xmin': lower_bound.x(),
                'xmax': upper_bound.x(),
                'ymin': lower_bound.y(),
                'ymax': upper_bound.y() 
}

# Check schema and table
if not schema or not table_name:
    raise GeoAlgorithmExecutionException('Must have a schema and a table name !')

if len(table_name) > 30:
    raise GeoAlgorithmExecutionException('Table name must be 30 characters max !')

f = open(str(output), 'w')
insert = u'-- Import '+layer.name()+u" in Oracle Spatial As " + table_name+u"...\n"
f.write(insert.encode('utf8'))

# Table creation:
if drop_table:
    insert = u'-- Drop the table:\n'
    insert = insert + u'DROP TABLE '+schema+"."+table_name + u';\n\n'
    f.write(insert.encode('utf8'))

insert = u"-- Create the table:\nCREATE TABLE "+schema+"."+table_name+ u'(\n'
f.write(insert.encode('utf8'))

# Field scanning to discover attributes length
field_length = {}

# Populating with the layer definition
for field in fields:
    if uppercase_field_names:
        name = field.name().upper()
    else:
        name = field.name()
    if field.type() == QVariant.String or field.type() == QVariant.Int:
        field_length[name] = str(field.length())
    elif field.type() == QVariant.Date or field.type() == QVariant.DateTime:
        field_length[name] = ''
    elif field.type() == QVariant.Double:
        field_length[name]=str(field.length())+","+str(field.precision())
    else:
        raise GeoAlgorithmExecutionException('Unknown Field type:'+field.typeName())

# Scanning features to determine true length
feats = processing.features(layer)
nFeat = len(feats)
for inFeat in feats:
    attrs = inFeat.attributes()
    i = -1
    for field in fields:
            i = i +1
            if uppercase_field_names:
                name = field.name().upper()
            else:
                name = field.name()
            if field.type() == QVariant.String:
                if (len(attrs[i])) > int(field_length[name]):
                    field_length[name] = str(len(attrs[i]))
            elif field.type() == QVariant.Int:
                if (len(str(attrs[i]))) > int(field_length[name]):
                    field_length[name] = str(len(str(attrs[i])))
            elif field.type() == QVariant.Double:
                current_len = int(field_length[name].split(',')[0])
                current_prec = int(field_length[name].split(',')[1])
                if (len(str(attrs[i]))) > current_len:
                    current_len = len(str(attrs[i]))
                if len(str(attrs[i]).split('.')) > 1:
                    if (len(str(attrs[i]).split('.')[1])) > current_prec:
                        current_prec = len(str(attrs[i]).split('.')[1])
                field_length[name] = str(current_len)+","+str(current_prec)

# Create the table
for field in fields:
    if uppercase_field_names:
        name = field.name().upper()
    else:
        name = field.name()
    insert = u"\t"+name
    if field.type() == QVariant.String:
        insert = insert + u" VARCHAR2("+field_length[name]+")"
    elif field.type() == QVariant.Date or field.type() == QVariant.DateTime:
        insert = insert + u" DATE"
    elif field.type() == QVariant.Int or field.type() == QVariant.Double:
        insert = insert +u" NUMBER("+field_length[name]+")"
    else:
        raise GeoAlgorithmExecutionException('Unknown Field type:'+str(field.type()))
    if field.name() == primary_key:
        insert = insert + u" CONSTRAINT PK_"+table_name[0:27]+u" PRIMARY KEY"
    insert = insert + u",\n"
    f.write(insert.encode('utf8'))

# Add Geometry definition
insert = u"\tGEOM MDSYS.SDO_GEOMETRY\n);\n"
f.write(insert.encode('utf8'))

# Remove Metadata
insert = u"\n-- Update Oracle Spatial Metadata\nDELETE FROM USER_SDO_GEOM_METADATA WHERE TABLE_NAME = '"+table_name+"' ;\nCOMMIT;\n\n"
f.write(insert.encode('utf8'))

# Add new metadata
insert = u"INSERT INTO USER_SDO_GEOM_METADATA (TABLE_NAME,COLUMN_NAME,DIMINFO,SRID)\n"
insert = insert + u"VALUES ('"+table_name+"','GEOM',\n"
insert = insert + u"\tMDSYS.SDO_DIM_ARRAY(\n"
insert = insert + u"\t\tMDSYS.SDO_DIM_ELEMENT('X',"+str(extent['xmin'])+u","+str(extent['xmax'])+u","+str(tolerance)+u"),\n"
insert = insert +u"\t\tMDSYS.SDO_DIM_ELEMENT('Y',"+str(extent['ymin'])+u","+str(extent['ymax'])+u","+str(tolerance)+u")),\n"
insert = insert + u"\t"+str(srid)+");\nCOMMIT;\n"
f.write(insert.encode('utf8'))

# Add table constraints
geom_types=[]
if geometry_type == 10:
    #Scan geometry to determine Geometry type
    feats = processing.features(layer)
    for inFeat in feats:
        the_type = inFeat.geometry().wkbType()
        if the_type == QGis.WKBPoint and 2001 not in geom_types:
            geom_types.append(2001)
        elif the_type == QGis.WKBMultiPoint and 2005 not in geom_types:
            geom_types.append(2005)
        elif the_type == QGis.WKBLineString and 2002 not in geom_types:
            geom_types.append(2002)
        elif the_type == QGis.WKBMultiLineString and 2006 not in geom_types:
            geom_types.append(2006)
        elif the_type == QGis.WKBPolygon and 2003 not in geom_types:
            geom_types.append(2003)
        elif the_type == QGis.WKBMultiPolygon and 2007 not in geom_types:
            geom_types.append(2006)

if geometry_type != 0:
    if geometry_type == 1:
        geom_types.append(2001)
    elif geometry_type == 2:
        geom_types.append(2005)
    elif geometry_type == 3:
        geom_types.append(2001)
        geom_types.append(2005)
    elif geometry_type == 4:
        geom_types.append(2002)
    elif geometry_type == 5:
        geom_types.append(2006)
    elif geometry_type == 6:
        geom_types.append(2002)
        geom_types.append(2006)
    elif geometry_type == 7:
        geom_types.append(2003)
    elif geometry_type == 8:
        geom_types.append(2007)
    elif geometry_type == 9:
        geom_types.append(2003)
        geom_types.append(2007)
    geom_types_str= [ str(sdf) for sdf in geom_types]
    geo_constraint=u"SDO_GTYPE = "+u" OR SDO_GTYPE = ".join(geom_types_str)
    insert = u"\n-- Geometry constraints\n"
    insert = insert + u"ALTER TABLE "+table_name+"\n"
    insert = insert+ u"\tADD CONSTRAINT GC_"+table_name[0:27]+" CHECK ("+geo_constraint+") ENABLE;\n"
    f.write(insert.encode('utf-8'))

# Add spatial_index:
if add_spatial_index:
    insert = u"\n-- Create Spatial Index\n"
    insert = insert + u"DROP INDEX GI_"+table_name[0:27]+";\n"
    insert = insert + u"CREATE INDEX GI_"+table_name[0:27]+" ON " +table_name+"(GEOM)\n"
    insert = insert +u"\tINDEXTYPE IS MDSYS.SPATIAL_INDEX\n"
    insert = insert + u"\tPARAMETERS (SDO_DML_BATCH_SIZE = 1');\n"
    insert = insert +u"COMMIT;\n"
    f.write(insert.encode('utf8'))


if inject_data:
    nElement = 0
    feats = processing.features(layer)

    for inFeat in feats:
        progress.setPercentage(int(100 * nElement / nFeat))
        nElement += 1
        # Starts INSERT
        insert = u"INSERT INTO "+table_name+ u" VALUES ("
        # Inject attributes
        attrs = inFeat.attributes()
        i = -1
        for field in fields:
            i = i +1
            if field.type() == QVariant.String:
                value = u"'"+attrs[i]+u"'"
            elif field.type() == QVariant.Date or field.type() == QVariant.DateTime:
                value = u"TO_DATE('"+attrs[i].toString('yyyy-MM-dd')+u"','YYYY-MM-DD')"
            else :
                value = str(attrs[i])
            insert = insert + value +u", "                       
        # Injects geometry in WKT format
        inGeom = inFeat.geometry().exportToWkt()
        insert = insert + u"SDO_GEOMETRY('" + inGeom + u"')"
        insert = insert + u");\n"
        f.write (insert.encode('utf8'))

# This is the end... 
f.close()

Conclusion

Le script n'est pas parfait mais il permet de générer du SQL conforme à Oracle Spatial et qui permettra à QGis de bien découvrir les couches nouvellement créées.

Quelques améliorations peuvent être imaginée. La première consiste à améliorer la gestion de l'import des données. En effet, même si j'ai indiqué les instructions pour réaliser un import des données et pas uniquement de la structure, cet import ne fonctionne pas vraiment dans la pratique. En effet, Oracle ne permet pas des requêtes SQL texte dont la taille est trop longue (limité à 4000 caractères je crois). Dans les faits, un import de polygones sera donc impossible alors qu'un import de points sera sans doute applicable. Il me faudrait trouver un autre moyen, notamment, profiter de ce que QGis offre pour créer une table et y injecter un contenu.

Enfin, on pourrait aller plus loin en permettant l'import direct dans Oracle Spatial, à l'instar du GéoAlgorithme qui existe pour PostGIS. Mais pour l'instant, il existe au moins un outil pour créer une structure de table géographique Oracle Spatial à partir d'une couche ouverte dans QGis. C'est déjà un bon début et ça permettra à ceux qui veulent lire le code de mieux comprendre comment on créé un GéoAlgorithme en Python...

Posted dim. 01 juin 2014 18:12:05

Voilà, c'est fait ! Suite à mon dernier article sur l'import du bâti du cadastre du Maine-et-Loire, j'avais eu pas mal de réactions qui m'enjoignaient de faire des corrections. En effet, mon import engendrait un certain nombre de problèmes dont le plus important était la remonté de nombreuses erreurs de croisement entre des voies de circulation et des bâtiments, ce qui, en temps normal, ne devrait pas se produire.

Cette erreur est essentiellement due au fait que le cadastre est généralement plus précis que le référentiel utilisé pour créer les routes dans OpenStreetMap. Pire, j'ai souvent vu des erreurs liées à une création de voie à partir des orthophotos de Bing. Très clairement dans certains coins, Bing a un vrai décalage par rapport au cadastre, non pas en précision mais en position pure: Bing est clairement à côté de la plaque dans ces lieux. Je pense sincèrement que le référentiel du cadastre est plus précis que Bing car il constitue un référentiel précis (au 1/5000).

Ces corrections ont clairement pris du temps mais pas autant que la réalisation de l'import pur, c'est déjà une bonne chose. La méthode que j'ai employée est assez simple. J'ai utilisé l'outil d'assurance qualité mis à disposition par la communauté OSM France et consultable sur http://osmose.openstreetmap.fr. Avec la requête suivante: il est possible d'avoir la liste des erreurs pour un utilisateur précis, selon le type d'erreur. Ensuite, j'ai utilisé la carte dynamique d'Osmose pour visualiser les problèmes à l'affichage. Dans un autre onglet, j'ouvrais une fenêtre OSM en mode modification avec l'éditeur iD.

L'intérêt d'avoir les deux fenêtres proches l'une de l'autre permet d'avoir en permanence la localisation géographique des problèmes à corriger tout en ayant la fenêtre de modifications sous la main. Je me suis donc appuyé uniquement sur le web pour réaliser toutes les corrections et j'ai trouvé ce mode de travail assez pratique. JOSM est plus complexe, l'import de données est beaucoup plus lent et ne m'a pas semblé pertinent pour réaliser un ensemble de corrections à la volée, réparties sur tout un département. J'ai donc pu constater à quel point l'éditeur iD est maintenant mûr et très fonctionnel. Grâce à lui, des modifications simples peuvent être réalisées, y compris si elles sont nombreuses, pour peu qu'on puisse s'astreindre à avoir un travail continu.

Je réserverai donc JOSM aux travaux d'import ou aux corrections massives sur un territoire raisonnable (une commune au maximum). De plus, je vous conseille également de faire attention au comportement de JOSM. Car ce qui m'a semblé étrange c'est que JOSM ne m'ait pas indiqué les erreurs de chevauchement. J'avais bien les erreurs liées à la géométrie (papillons ou bâtiments se chevauchant par exemple) mais pas les problèmes d'intersection avec les voies. Pourtant JOSM sait gérer ce genre de problèmes. Si vous réalisez des imports de bâti, je vous conseille sincèrement de trouver le moyen de contrôler ces erreurs sous peine de devoir faire comme moi et de retravailler à posteriori.

J'ai été assez surpris par le fait qu'il y avait beaucoup de monde pour me faire remarquer que mon import posait problème mais finalement assez peu de personnes pour faire les corrections ou pour me proposer leur aide. En effet, le petit millier d'erreurs est resté assez volumineux jusqu'à ce que je m'y attelle. Je dois toutefois remercier ceux qui ont contribué et que je n'ai peut-être pas remarqué. Dans tous les cas, mon travail d'import est terminé et c'est une bonne chose...

Posted mar. 17 juin 2014 20:12:00 Tags: