QGIS Processing Grass Provider explanations
Last year (2017) I worked on the QGIS Grass Processing Provider. With
time, I finally found how it worked. It is not so complicated but it
is not so easy to understand because there is a "trick" on what is
called the ext mechanism. Sometimes, developpers contact me by
asking for help on this mechanism.
Here is my last email response about this subject. I also put it on my blog for people who want to understand or work on this part of QGIS.
Hello,
it's been a little bit more than one year since I have "retired" from
the QGIS project. It's just because I do have other projects that I am
working for and I had to cut down from some projects.
Well, I can only tell you about the GRASS Processing from QGIS 2.18
(but I bet that this is still the same in QGIS 3).
If you want to use the QGIS GRASS plugin, it is completely different
and I can't help you on this part.
All you have to know is that:
* The GRASS provider for QGIS Processing is just a graphical interface
in QGIS towards GRASS.
* There is an on-the-fly creation of a GRASS database in a temporary
directory. Then all the input layers/rasters are imported (or
referenced) into this GRASS database.
* When you launch an algorithm from the GRASS provider, QGIS launch
(via the shell or a system call) GRASS command line with the right
parameters onto the temporary GRASS database. You can read this
command line in the log.
* If the algorithm returns some results, those are extracted from the
GRASS temporary database and converted into layers/rasters (and
directly opened into QGIS if option has been enabled). Then the
GRASS temporary database is deleted.
To add your add-on in the GRASS provider, you need to:
* Build a description file in
python/plugins/processing/algs/grass7/description/. Those files just
describe the different parameters of your algorithm (add-on). You
can read some of those files to understand the different
options. The description file filename must respect the algorithm
name.
* If your algorithm doesn't require something else than: import input
data/launch the algorithm/extract data, then you are done. There is
nothing more to do than the description file.
* If your algorithm need more actions like a preliminary check on
data, more than one GRASS command involved or a transformation on
output data before extracting them, then you need to use the ext
mechanism.
The 'ext' mechanism is a (some say kind of dirty) way to add per
algorithm additional logic. There are 4 different levels where you can
add logic:
- Checking the input parameters, if you want to verify that two
excluding options have been both enabled.
- Processing inputs import: if you need to do more than importing
input layers.
- Processing the command by itself: if you need to chain more than one
GRASS command for your algorithm.
- Processing the outputs export: if you need to do special things
before exporting layers or if you need special export methods.
Whenever you want to add some logic on one (or more) level(s), you
have to create a .py file named accordingly to the algorithm name
(modify '.' with '_'), into
python/plugins/processing/algs/grass7/ext. Then you need to create
methods that will respect names:
- Input parameters: checkParameterValuesBeforeExecuting
- Inputs import: processInputs
- Command: processCommand
- Outputs: processOutputs
If there is a Python file with the algorithm name in the ext
directory, methods will be imported from the file and run instead of
the common methods (there are "standard"
processCommand/processInputs/processOutputs/checkParameterValuesBeforeExecuting
in the code of the GRASS provider for QGIS Processing, in
python/plugins/processing/algs/grass7/Grass7Algorithm.py).
If we take the example of v.what.rast, there is an ext file:
ext/v_what_rast.py. In this file there is a processCommand method. It
just launch the standard processCommand but with the
delOutputs option set to True (we do not want to have standard
outputs). Then there is also a customized processOutputs which export
the input vector as an output for QGIS. We need to do this because
v.what.rast modify values directly into the input vector layer. It
doesn't generate dedicated output so we have to build this output
ourself.
Clearly, if you want to do special things in the ext mechanism, you
will need to read (and understand) the GRASS provider code standard
methods in Grass7Algorithm.py.
I hope that these explanations will ease your work...
Tracer des isochrones avec des données et des logiciels libres
Introduction
Ça fait bien longtemps que je n'ai pas discuté de SIG sur ce blog. Il faut dire que j'en fait un peu tous les jours, dans le cadre de mon activité professionnelle. Mais de temps en temps, c'est bien de se remettre dans un cadre plus libre.
Cet article a pour objectif de présenter une méthode pour tracer des isochrones: des courbes qui indiquent à quelle distance on peut aller dans une enveloppe de temps à partir d'un ou plusieurs points de départ différents.
Nous allons utiliser des outils libres mais également des données libres, avec toutes leurs contraintes. Notre objectif est de voir , en partant des gares d'un département, quelles sont les zones où l'on peut se rendre en voiture pendant un laps de temps donné. Nous allons retenir les plages de 5 minutes/10 minutes/18 minutes et 20 minutes.
Cet article n'a pas pour but d'être une référence de ce genre de calcul. Il présente juste une méthode parmi d'autres. Elle peut d'ailleurs se révéler très approximative et non adaptée à certains besoins. Voilà vous êtes prévenus...
Capture de données géographiques
Qu'est-ce-qu'il nous faut ?
Avant de procéder, il faut disposer de données géographiques. C'est la base. Voici une petite liste de ce dont nous avons besoin:
- Le tracé des routes, c'est la base. C'est ce qui va constituer notre réseau de déplacement.
- Les limites de comunes et de département pour pouvoir se situer plus facilement.
- Les zones résidentielles pour compléter les indications de vitesse des routes.
- Les gares et arrêts SNCF, ce seront nos points de départ ou d'arrivée. C'est ce qui va constituer les noeuds d'intérêt de notre réseau de déplacement.
Pour nos besoins, nous allons utiliser la base de données OpenStreetMap. Cette dernière est libre de droits et elle est assez complète même si, comme nous le verrons, elle pourrait être améliorée fortement.
Les données via l'API overpass
L'intérêt d'OpenStreetMap, c'est qu'on peut effectuer des requêtes via une API. Le language est assez complexe mais on peut vraiment récupérer ce qu'on veut.
Un moyen simple d'effectuer des requêtes est d'utiliser le requêteur overpass-turbo. Il permet de faire des requêtes à partir d'un navigateur web et, surtout, de visualiser les résultats avant de pouvoir les télécharger.
Voici un exemple commenté de requête pour récupérer des routes:
// On veut une sortie au format JSON
// Au bout de 155 secondes, on arrête la requête.
[out:json][timeout:155];
// On indique que la zone de recherche est identifié par le polygone
// dont le nom est "Le Loroux-Bottereau", une commune du département
// de Loire-Atlantique.
{{geocodeArea:"Le Loroux-Bottereau"}}->.searchArea;
(
// On indique qu'on veut récupérer les routes
way["highway"](area.searchArea);
);
// Les options pour la sortie.
out body;
>;
out skel qt;
Avec cette requête, on peut récupérer les routes d'une commune. Néanmoins, pour notre exemple, il nous faut télécharger les routes de toutes les communes d'un département. Le faire à la main est forcément fastidieux.
Certains penseront à modifier le nom de recherche en y placant celui du département. Néanmoins, ils atteindront les limites de l'API. En effet, les serveurs qui fournissent le service de l'API sont partagés et ne peuvent passer trop de temps à effectuer une requête.
Dépasser les limites de l'API
Pour récupérer l'ensemble de nos routes, nous avons deux solutions:
- Télécharger une archive toute faite.
- Scripter la récupération des communes par un script shell.
Et si on faisait quand même un petit script ?
Même si c'est la méthode la moins rapide, ça vaut toujours le coup de disposer d'une méthode alternative, par exemple pour le cas où vos données ne sont pas disponibles dans le téléchargement en masse.
Concrètement, l'avantage de l'API overpass, c'est que c'est une API
HTTP. On peut donc l'utiliser avec un simple client HTTP en ligne de
commande comme wget ou curl (ce dernier a ma préférence).
Dans notre cas, nous avons deux choses à régler: * Déterminer le code de la zone de recherche: sur le site overpass-turbo, on peut utiliser un nom comme zone de recherche mais, en vérité, une première requête est lancée sur nominatim pour transformer ce nom en identifiant de zone. * Lancer une boucle de requêtes qui sont suffisamment lentes pour ne pas effrayer les serveurs de l'API.
Voici un exemple de script shell commenté qui répond à ces problèmes:
#!/bin/bash
echo "Récupération des routes des communes..."
# Notre fichier d'entrée est un fichier csv qui contient deux
# colonnes: le code osm_id et le nom de la commune.
# Vous pouvez créer ce fichier en téléchargeant les communes du
# département via l'API overpass et en utilisant QGIS pour fabriquer le CSV.
INPUT=communes.csv
OLDIFS=$IFS
IFS=,
[ ! -f $INPUT ] && { echo "$INPUT file not found"; exit 99; }
# Le début de notre boucle
while read comid comname
do
# On complémente à 8 chiffres l'identifiant OSM
comid=$(printf "%08d" $comid)
echo "Dowloading roads for $comname... ($comid)"
# On précise notre requête OverPass
request="[out:xml][timeout:155];area(36"${comid}')->.searchArea;(way["highway"](area.searchArea););(._;>;);out
body;'
# et on la joue
curl 'https://overpass-api.de/api/interpreter' -H 'Content-Type: application/x-www-form-urlencoded; charset=UTF-8' -H 'Origin:https://overpass-turbo.eu' --data "$request" -o "roads_${comname}.osm"
# Une méthode bourrine pour gérer les limite de l'API: on attend 10 secondes entre chaque requête.
sleep 10
done < $INPUT
IFS=$OLDIFS
# Un script shell se termine toujours par une valeur de retour !
exit 0
Limites de département
Voici la requête overpass qui permet de récupérer la limite du départment:
[out:xml][timeout:155];
{{geocodeArea:"Loire-Atlantique"}}->.searchArea;
(
rel["admin_level"="6"](area.searchArea);
);
(._;>;);out body;
Limites des communes
Voici la requête overpass qui permet de récupérer les limites de communes:
[out:xml][timeout:155];
{{geocodeArea:"Loire-Atlantique"}}->.searchArea;
(
rel["admin_level"="8"](area.searchArea);
);
(._;>;);out body;
Pour l'extraction depuis OSM, utiliser grabCommunes qui utilise l'API Overpass. Pour l'import
Les gares
Voici la requête overpass qui permet de récupérer les gares et arrêts de train du départment:
[out:json][timeout:155];
{{geocodeArea:"Loire-Atlantique"}}->.searchArea;
(
node["railway"="station"](area.searchArea);
node["railway"="stop"](area.searchArea);
);
(._;>;);out body;
Les zones résidentielles
Voici la requête overpass qui permet de récupérer les zones résidentielles:
[out:json][timeout:155];
{{geocodeArea:"Loire-Atlantique"}}->.searchArea;
(
way["landuse"="residential"](area.searchArea);
);
(._;>;);out body;
Traitement des données
Chargement
Avant de pouvoir réaliser des traitements et des analyses sur les données que nous venons de récupérer, il faut pouvoir les stocker dans un endroit adapté. En effet, nous disposons de fichiers XML qui sont certes utilisables directement dans des logiciels comme QGIS, mais qui se révèlent peu performants (pour cause d'absence d'index spatial par exemple).
De plus, nous allons devoir réaliser des traitements spécifiques qui ne sont pas accessibles à tous les formats de stockage. Pour pouvoir aller plus vite, autant prendre ce qui est le plus performant, le plus complet: PostGIS.
Je vous laisse le soin de monter un serveur PostgreSQL avec
l'extension PostGIS sur votre bécane. Pour ma part, je travaille avec
un serveur qui écoute sur localhost, la base de données se nomme
geobase, l'utilisateur se nomme geoadmin et il n'a pas besoin de
mot de passe pour se connecter (c'est mal mais rapide).
Pour le chargement des XML, j'utilise ogr2ogr sur cette forme:
ogr2ogr -overwrite -f "PostgreSQL" "PG:host=localhost dbname=geobase user=geoadmin" -t_srs EPSG:2154 -lco "OVERWRITE=YES" -lco "SPATIAL_INDEX=YES" -lco "GEOMETRY_NAME=geom" -nlt MULTILINESTRING -nln communes -progress initCommunes.osm lines
Cette commande permet de créer une couche PostGIS nommée communes dans
la base indiquée, dans la projection Lambert93 (2154), avec des index
spatiaux, où le champ de géométrie se nomme geom.
Si vous avez plusieurs fichiers à charger, je vous conseille d'utiliser la technique suivante: * Renommer un des fichiers que vous avez à importer en initCommunes.osm * Utilisez le script suivant:
#!/bin/bash
# Import initial, permet de créer la structure
ogr2ogr -overwrite -f "PostgreSQL" "PG:host=localhost dbname=geobase user=geoadmin" -t_srs EPSG:2154 -lco "OVERWRITE=YES" -lco "SPATIAL_INDEX=YES" -lco "GEOMETRY_NAME=geom" -nlt MULTIPOLYGONS -nln communes -progress initCommunes.osm multipolygons
# Boucle d'ajout !
find ./ -name 'communes_*.osm' | while read com
do
echo "$com"
ogr2ogr -append -f "PostgreSQL" "PG:host=localhost dbname=geobase user=geoadmin" -t_srs EPSG:2154 -nlt MULTIPOLYGONS -nln communes -progress "$com" multipolygons
done
exit 0
Attention, pour les données des gares, je vous conseille de supprimer les données en double: certaines gares ont également des points d'arrêt.
Amélioration des données de routes: phase 1
Les données brutes sont chargées mais il reste néanmoins à les retravailler:
- D'abord, il y a des chances que les données de routes soient en doublons dans certains cas.
- Ensuite, il faut faire un peu de ménage dans les champs de cette couche.
- Il faut croiser cette donnée avec les zones résidentielles.
- Enfin, il faut s'assurer de disposer d'un réseau topologique correct.
Pour ces opérations de corrections, nous allons utiliser PostGIS et GRASS.
Dans un premier temps, il faut faire un peu de ménage, dans PostGIS:
-- Nettoyage des données de routes
---- Suppression des doublons géographiques
DELETE FROM routes WHERE ogc_fid IN (select ogc_fid from (
SELECT ogc_fid, ROW_NUMBER() OVER(PARTITION BY geom ORDER BY ogc_fid asc) AS Row,
geom FROM ONLY public.routes
) dups where dups.Row > 1
order by ogc_fid);
---- Un peu de ménage, on supprime les routes non utiles dans notre étude
DELETE FROM routes WHERE highway IN
('track', 'bus_stop', 'bus_guideway', 'construction', 'crossing', 'cycleway',
'disused', 'footway', 'proposed', 'raceway', 'rest_area', 'service', 'access',
'bridleway', 'path', 'pedestrian', 'platform','steps', 'road');
---- Modification du champ nom s'il est vide
UPDATE routes SET name = substring(other_tags from '"ref"=>"([^\"]+)".*$')
WHERE name IS NULL
AND other_tags ~ '.*"ref"=>".*';
---- Un peu de ménage dans les champs
ALTER TABLE routes DROP COLUMN IF EXISTS waterway;;
ALTER TABLE routes DROP COLUMN IF EXISTS aerialway;
ALTER TABLE routes DROP COLUMN IF EXISTS z_order;
ALTER TABLE routes DROP COLUMN IF EXISTS barrier;
ALTER TABLE routes DROP COLUMN IF EXISTS man_made;
ALTER TABLE routes ADD COLUMN maxspeed NUMERIC(10,2);
---- Pré-affectation des vitesses en fonction du type de route
---- Récupération de la donnée OSM
UPDATE routes SET maxspeed = substring(other_tags from '"maxspeed"=>"([^\"]+)".*$')::NUMERIC(10,2)
WHERE other_tags ~ '.*"maxspeed"=>".*'
AND substring(other_tags from '"maxspeed"=>"([^\"]+)".*$') ~ '^\d+$';
---- Les cas spéciaux sont à 50 km/h
UPDATE routes SET maxspeed = 50.0
WHERE other_tags ~ '.*"maxspeed"=>".*'
AND maxspeed IS NULL;
UPDATE routes SET maxspeed = 130.0 WHERE highway IN ('motorway') AND maxspeed IS NULL;
UPDATE routes SET maxspeed = 90.0 WHERE highway IN ('motorway_link', 'trunk_link') AND maxspeed IS NULL;
UPDATE routes SET maxspeed = 50.0 WHERE highway IN ('residential') AND maxspeed IS NULL;
-- Croisement avec les zones résidentielles
---- Aggrégation des zones résidentielles
DROP TABLE IF EXISTS agg_resi;
CREATE TABLE agg_resi AS SELECT 1 as gid, ST_Union(p.geom) as geom from residentials p;
CREATE INDEX agg_resi_geomidx
ON public.agg_resi
USING gist
(geom);
-- On créé une colonne pour indiquer qu'on veut supprimer des lignes
ALTER TABLE routes DROP COLUMN IF EXISTS todelete;
ALTER TABLE routes ADD COLUMN todelete BOOLEAN DEFAULT False;
-- On remplit cette colonne avec l'information
UPDATE routes a SET todelete = True
FROM residentials p
WHERE ST_Intersects(a.geom, p.geom)
AND a.highway IN ('primary','secondary','tertiary', 'unclassified');
-- On ajoute les routes découpées
--- Attention, cette requête prend beaucoup de temps...
INSERT INTO routes (osm_id, name, highway, other_tags, maxspeed, todelete, geom)
SELECT a.osm_id, a.name, a.highway, a.other_tags, a.maxspeed, a.todelete,
ST_Multi(a.clipped_geom)
FROM (
SELECT l.osm_id, l.name, l.highway, l.other_tags, 50::integer AS maxspeed, False As todelete,
(ST_Dump(ST_Intersection(p.geom, l.geom))).geom As clipped_geom
FROM routes l
INNER JOIN agg_resi p ON ST_Intersects(l.geom, p.geom)
WHERE l.highway IN ('primary','secondary','tertiary', 'unclassified')
UNION
SELECT l.osm_id, l.name, l.highway, l.other_tags, 80::integer AS maxspeed, False As todelete,
(ST_Dump(ST_Difference(l.geom, p.geom))).geom As clipped_geom
FROM routes l
INNER JOIN agg_resi p ON ST_Intersects(l.geom, p.geom)
WHERE l.highway IN ('primary','secondary','tertiary', 'unclassified')
) a
WHERE ST_Dimension(a.clipped_geom) = 1;
-- On retire les routes à supprimer
DELETE FROM routes WHERE todelete = True;
-- On supprime la colonne todelete
ALTER TABLE routes DROP COLUMN todelete;
-- On vire la table d'aggrégation
DROP TABLE IF EXISTS agg_resi;
Amélioration des données de routes: phase 2
Ensuite, il nous faut créer un vrai réseau topologique. En effet, ce qui sort d'OSM est un ensemble de routes qui ne se coupent pas forcément au niveau des intersections. Globalement, ça peut poser des problèmes pour les traitements futurs de GRASS.
Pour y pallier, il va nous falloir faire un peu de ménage avec GRASS. Vous devez donc disposer d'une DB GRASS et avoir GRASS7 de lancé. Nous allons juste lancer un script de nettoyage qui a la forme suivante:
#!/bin/bash
# Nettoyage du réseau routier dans GRASS
## On effectue un lien vers les données PostGIS
echo "Import des données de route depuis PostGIS..."
g.remove -f type=vector name=routes
g.remove -f type=vector name=c_routes
v.in.ogr --overwrite -o input="PG:dbname=geobase user=geoadmin host=localhost" layer=routes output=routes type=line
# Nettoyage des routes (on recalcule la topologie et on vire les doublons)
v.clean --overwrite -c input=routes output=c_routes type=line tool=break,rmline,rmdupl
# On exporte en GPKG car GRASS a un problème pour écrire directement
# dans la base PostgreSQL
v.out.ogr --overwrite input=c_routes output="data.gpkg" format=GPKG layer=1 type=line output_layer=c_routes
# puis on utilise ogr2ogr pour faire le transfert
ogr2ogr -overwrite -f "PostgreSQL" "PG:host=localhost dbname=geobase user=geoadmin" -t_srs EPSG:2154 -lco "OVERWRITE=YES" -lco "SPATIAL_INDEX=YES" -lco "GEOMETRY_NAME=geom" -nlt MULTILINESTRING -nln croutes -progress data.gpkg c_routes
exit 0
En sortie, on dispose de la couche croutes dans PostGIS. Il nous reste à calculer les coûts de déplacement.
Amélioration des données de routes: phase 3
Maintenant qu'on dispose d'un réseau topologique correct, il reste une dernière étape de calcul des coûts de déplacement. Comme nous souhaitons réaliser des isochrones, le coût va être exprimé dans une unité de temps (secondes dans notre cas).
J'ai pris l'hypothèse fortement discutable de prendre un temps de trajet par tronçon égal au temps de parcours à la vitesse maximale du tronçon.
Néanmoins, pour pondérer un peu les résultats, j'ai délibérément abaissé les vitesses en suivant les règles suivantes:
- Vitesse limitée à 90 km/h => 80 km/h
- Vitesse limitée à 80 km/h => 70 km/h
- Vitesse limitée à 50 km/h => 40 km/h
- Vitesse limitée à 30 km/h => 25 km/h
Mettons tout ça en musique via PostGIS:
-- Calcul des coûts de déplacement
-- On va dire que toutes les routes qui n'ont pas de maxspeed sont à 50 km/h pour les non classées
UPDATE croutes SET maxspeed = 50.0 WHERE maxspeed IS NULL AND highway IN ('unclassified', 'living_street');
-- On va dire que toutes les routes qui n'ont pas de maxspeed sont à 80 km/h pour les non classées
UPDATE croutes SET maxspeed = 80.0 WHERE maxspeed IS NULL AND highway NOT IN ('unclassified', 'living_street');
-- On abaisse les vitesses suivant les règles sus-citées
UPDATE croutes SET maxspeed = 80.0 WHERE maxspeed = 90.0;
UPDATE croutes SET maxspeed = 70.0 WHERE maxspeed = 80.0;
UPDATE croutes SET maxspeed = 50.0 WHERE maxspeed = 50.0;
UPDATE croutes SET maxspeed = 25.0 WHERE maxspeed = 30.0;
-- Enfin, on calcule les coûts de déplacement.
ALTER TABLE croutes DROP COLUMN IF EXISTS cost;
ALTER TABLE croutes ADD COLUMN cost NUMERIC(10,2);
UPDATE croutes set cost = (ST_Length(geom)/(maxspeed/3.60)) where maxspeed > 0;
Ok, nous avons de quoi calculer des isochrones.
Calcul des isochrones
Ce calcul va se faire dans GRASS, non pas parce que c'est l'optimum, mais juste parce que c'est la seule méthode que je connaisse, en plus de PgRouting.
Pour le calcul des isochrones, nous allons utiliser nos données de routes qui vont constituer le réseau et nos données de gares qui vont constituer des centres de départ sur ce réseau.
GRASS dispose de fonctions dédiées au réseau (les fonctions v.net), que nous allons utiliser.
Le script suivant contient ce qu'il faut pour réaliser le calcul d'isochrones. Il doit être lancé dans l'environnement shell de GRASS:
#!/bin/bash
# Calcul d'isochrones dans GRASS
# On relie les données utiles dans la base GRASS sans faire d'import.
# C'est inutile dans notre cas car les fonctions de réseau vont dupliquer les données à leur sauce.
echo "Import des données depuis PostGIS..."
v.external -o --quiet --overwrite input="PG:host=localhost user=geoadmin dbname=geobase" layer=gares output=gares
v.external -o --quiet --overwrite input="PG:host=localhost user=geoadmin dbname=geobase" layer=croutes output=croutes
# Nettoyage des routes
# v.clean --overwrite -c input=routes output=c_routes type=line tool=break,rmline,rmdupl
# Construction du réseau
## Ici, on va construire un réseau au sens GRASS du terme. On part de notre couche de routes, on la convertie en réseau et on y adjoint les centres (les gares). Elles sont ajoutées au réseau par un trait de base.
echo "Construction du réseau..."
v.net --quiet -s --overwrite input=croutes points=gares output=thenetwork operation=connect arc_type=line threshold=50
# lier les bases de données des routes et des noeuds (uniquement si v.in.ogr):
#v.db.connect -o map=thenetwork@PERMANENT table=gares layer=2
# v.db.connect -o map=network@PERMANENT table=routes layer=1
# Calcul des isochrones
echo "Calcul des isochrones..."
## Ici, on prend toutes les gares et on calcule les isochrones sur des limites de coûts de 300, puis de 600 puis de 1080 puis de 1200.
## Cette opération casse le réseau selon les limites. Comme notre coût est un temps, 300 correspond donc à 5 minutes (300 secondes).
v.net.iso -u --overwrite input=thenetwork@PERMANENT output=isochrones@PERMANENT center_cats=1-1000 costs=300,600,1080,1200 arc_column=cost
# Nettoyage des valeurs à -1 dans les isochrones
echo "Suppression des valeurs inutiles..."
v.edit --quiet map=isochrones tool=delete layer=1 type=line where="center=-1 or isonr>=5"
# Exporter les données dans PostGIS (rapide) mais doit être défini dans db.connect.
echo "Export des données..."
v.out.ogr --overwrite -c input=isochrones output="PG:dbname=geobase" format=PostgreSQL layer=1 type=line lco="OVERWRITE=YES,SCHEMA=public,GEOMETRY_NAME=geom,FID=fid" output_layer=isochrones
exit 0
Créer les polygones des zones
Vous devriez pouvoir afficher les polylignes des isochrones. Concrètement, chaque coût cumulé a été injecté dans chaque tronçon de route suivant les 4 classes énoncées plus haut (300,600,1080,1200).
Il reste à générer des polygones à partir de ces éléments. Pour ça, le plus simple consiste à utiliser PostGIS avec les requêtes qui suivent:
-- Calcul des zones d'isochrones
---- On commence par créer la tables des zones iso
DROP TABLE IF EXISTS isozones;
CREATE TABLE public.isozones
(
fid serial,
center integer,
classe integer,
geom geometry(MultiPolygon, 2154),
CONSTRAINT isonzones_pk PRIMARY KEY (fid)
)
WITH (
OIDS=FALSE
);
ALTER TABLE public.isozones
OWNER TO geoadmin;
CREATE INDEX isozones_geom_idx
ON public.isozones
USING gist
(geom);
-- Remplissage de la table des zones
INSERT INTO isozones (center, classe, geom) SELECT center, isonr, ST_Multi(ST_ConvexHull(ST_Collect(geom))) as geom FROM isochrones GROUP BY center, isonr;
-- Zones fusionnées
---- Ensuite, on va aggréger les zones par classe
DROP TABLE IF EXISTS merged_isozones;
CREATE TABLE public.merged_isozones
(
fid serial,
classe integer,
geom geometry(MultiPolygon, 2154),
CONSTRAINT merged_isozones_pk PRIMARY KEY (fid)
)
WITH (
OIDS=FALSE
);
ALTER TABLE public.merged_isozones
OWNER TO geoadmin;
CREATE INDEX merged_isozones_geom_idx
ON public.merged_isozones
USING gist
(geom);
-- Remplissage de la table des zones
TRUNCATE TABLE merged_isozones;
INSERT INTO merged_isozones (classe, geom) SELECT classe, ST_Multi(ST_Union(geom)) as geom FROM isozones GROUP BY classe;
-- Découpage manuel des classes
---- Je sais, c'est crade mais c'est rapide.
UPDATE merged_isozones SET geom = ST_Multi(ST_Difference(b.geom, a.geom))
FROM (SELECT geom FROM merged_isozones WHERE classe = 1) a,
(SELECT geom FROM merged_isozones WHERE classe = 2) b
WHERE classe = 2;
UPDATE merged_isozones SET geom = ST_Multi(ST_Difference(b.geom, a.geom))
FROM (SELECT ST_Union(geom) as geom FROM merged_isozones WHERE classe IN (1,2)) a,
(SELECT geom FROM merged_isozones WHERE classe = 3) b
WHERE classe = 3;
UPDATE merged_isozones SET geom = ST_Multi(ST_Difference(b.geom, a.geom))
FROM (SELECT ST_Union(geom) as geom FROM merged_isozones WHERE classe IN (1,2,3)) a,
(SELECT geom FROM merged_isozones WHERE classe = 4) b
WHERE classe = 4;
Affichage dans QGIS
Les données sont maintenant disponibles pour visualisation dans QGIS. En guise de fond de plan, vous pouvez utiliser les tuiles TMS d'OpenStreetMap:
https://a.tile.openstreetmap.org/{z}/{x}/{y}.png
Avec les éléments présentés précédemment, j'arrive à constituer la carte suivante:
Tout ça pour ça ? Oui mais maintenant vous savez faire. A vous de l'adapter à vos besoins.
Conclusions
Bon, c'était assez compliqué et ce, pour plusieurs raisons. La première c'est que les données d'OpenStreetMap ne sont pas encore assez fines, notamment au niveau des limites de vitesse. J'ai dû batailler pour trouver quelques techniques pour les déduires d'autres couches.
Par ailleurs, les routes OpenStreetMap ne forment pas du tout un réseau topologique (enfin du moins, en sortie de l'API overpass) ce qui pose de nombreux problèmes pour GRASS. Nous arrivons à contourner ces éléments mais cela impose d'ajouter quelques requêtes.
Enfin, ça fait longtemps que je n'ai pas vraiment fait de SIG libre: je suis donc un peu rouillé et il doit y avoir des techniques pour aller plus vite, avec d'autres opérateurs.
Le choix de PostGIS peut sembler bizarre de prime abord, essentiellement parce que finalement, tous les traitements se font sur la même machine. Au départ, j'avais commencé avec le format standard GPKG qui est vraiment bien géré par QGIS et GDAL. Mais finalement, vu les traitements imposés pour le nettoyage des données, je me suis retrouvé à court de fonctions géométriques dans Spatialite (GPKG est une surcouche de Spatialite). Je me suis rabattu vers PostgreSQL/PostGIS qui a effectivement tenu ses promesses.
Au final, on obtient des isochrones avec tous les détails pour rejouer les calculs avec d'autres classes de coût de déplacement. On pourrait aussi aller plus loin en recalculant les coûts en fonction de la présence d'éléments de signalisation qui imposent ou non un arrêt.
Mais d'ici là, et en restant approximatif, on obtient des zones pas si délirantes que ça en comparant avec ce que peuvent produire OSRM ou même Google Maps dans leur calcul de distance.
Pour terminer, on aurait pu s'affranchir de GRASS si nous avions utilisé PgRouting. Mais je souhaitais aussi montrer que le SIG GRASS avait encore des atouts sur la partie des réseaux.
Introduction
This article is the last of my series of articles about QGIS forms. You can read the previous one here.
Today, we will focus on n,n relations. n,n is a database syntax to tell that one object of a table can have multiples values linked to another reference table. On a database schema, it looks like the following:
You can see that n,n relations are materialized by two 1,n relations, involving three tables:
- The first data table (ANALYSIS).
- The second data table (PESTICIDE).
- A relation table (ANALYSIS_PESTICIDE) that make the connexion between the two tables mentionned above.
With such a mechanism, we are able to link multiple pesticides to multiple analysis without having to store a list into the ANALYSIS table. Instead, the central table (ANALYSIS_PESTICIDE) is used to make a link. Everytime you would like to implement n,n relations, just think about this intermediate table.
For the moment, QGIS doesn't support n,n tables on forms (but it supports 1,n with sub-forms). There is no control dedicated to that. But we can code it !
User Interface
What to do ?
Before diving into code, we need to solve the UI problem. What control are we going to use ? As we are not constrained by QGIS existing controls, we can imagine the following:
- The main source of data is ANALYSIS: we want to add multiple pesticides from one analysis.
- We only need the relevant information: which pesticides have been measured during the analysis.
- If you have a catalog of ten thousand of pesticides, there is no need to show the whole list on the analysis form.
- So we need a dedicated dialog for selecting pesticides. In this sub-form, the list of all the pesticides will be presented to the user in order to make a choice.
- The pesticides table will store all informations about pesticides. We need to have a list of the pesticide's names and be able to choose multiple pesticides entries.
- The simplest way to achieve this is to use a QListWidget with a checkbox for each pesticide entry. If the checkbox is checked, this analysis has this pesticide.
- We need a "search engine" to quickly find pesticides from their names if the PESTICIDE table is huge.
- Once selected, we need to only show the selected pesticides into analysis form.
Main form: analysis form
Here is a mockup of the analysis form:
You can see that the last form control is a bit special: it is our QListWidget which lists all the pesticides that have been found in the displayed analysis. The list only shows the relevant information and when the control is not long enough, you have a vertical scrollbar. Elements of the list are selectable (multi or mono depending on QListWidget attributes) and you can copy/paste them in the clipboard. Whenever you need to have some information about the analysis, everything is displayed on only one form.
This form can't be autogenerated by QGIS because the QListWidget is mandatory and there is no field to hold pesticides values in ANALYSIS table. For n,n relations you have to use a custom ui form with qt4-designer like the following:
Furthermore, printing information in the QListWidget of this form needs some code for a dedicated function in Python just to retrieve the good results from the ANALYSIS_PESTICIDE table. We will study this in the code part below.
What about editing pesticides into this analysis ?
Pesticide form
When you click on the "Modify" button, the following dialog will be printed:
The dialog box is very simple: on the top, you have a QLineEdit which will be used to type the name of the pesticide you want. The main control of the dialog is a QListWidget with the name of all the pesticides. There is a checkbox to add or remove pesticides to the analysis. Checking a box will add the pesticide into the ANALYSIS_PESTICIDE table, unchecking will delete it from the table. With this dialog, you can add or delete as many pesticides you want for one analysis without bothering with the other controls.
Whenever your modifications are done, results need to affect the ANALYSIS_PESTICIDE table and it also needs some dedicated code.
Using QGIS relations and value relations
Now that the concepts of the UI part have been elaborated, we need to go further. Our approach seems to be good with one n,n relation. But imagine that you try to build a true complex GIS application that involves about 50 n,n relations. you can't put everything into code, it will take too much time to develop and to maintain. You will also make a lot of mistakes to try to keep the names of all the controls into Python code. So we need to be a little bit more generic.
QGIS has already a mechanism to handle 1,n relations: it is called "Relations". Relations are a way for QGIS to know that a table is linked to another. It is used to show sub-forms inside a parent form. So, could we try to use two 1,n relations and deal with it into code ? I have tried this but there is something more efficient. Creating two relations (first from ANALYSIS to ANALYSIS_PESTICIDE and second from PESTICIDE TO ANALYSIS_PESTICIDE) seems to be the good way but you have to remember what we want. We would like to display a list of pesticides names in the control and store the ID_PESTICIDE into ANALYSIS_PESTICIDE and there is nothing in a QGIS relation to tell that you want to display a field.
But there is something inside QGIS to deal with and it's called "Value Relation". When you define a Value Relation for a field, you are linking values from another layer and you can choose what field to show and what field is the ID. So instead of creating two relations, we could do the following:
- Create only one relation: from ANALYSIS to ANALYSIS_PESTICIDE.
- Create a Value Relation control for ID_PESTICIDE of ANALYSIS_PESTICIDE towards ID_PESTICIDE of table PESTICIDE and show PESTICIDE.NAME (this is the field we want to display in the n,n sub dialog).
Here is the definition of the relation:
The name/id of the relation should be the same than the QListWidget of the custom .ui file.
Here is the definition of the Value Relation:
This is a classic Value Relation configuration. It is made in the ANALYSIS_PESTICIDE table on the ID_PESTICIDE attribute.
The relation is used as the following:
- the code launched when the ANALYSIS form is opened knows that the layer of the form is ANALYSIS.
- with this name, it is easy to search inside the project relations (there is an API for that) where ANALYSIS is the parent layer of another.
- code will take the name of the child layer (ANALYSIS_PESTICIDE) and extract the shared attributes (ID_ANALYSIS).
- Now, we know from what table we must read the results to update the form (this will be ANALYSIS_PESTICIDE) and we also know what is the attribute to filter (ID_ANALYSIS) to have the corresponding value of the analysis that is displayed in the form.
The Value Relation is used as the following:
- We already know that the intermediate table is ANALYSIS_PESTICIDE.
- We search for its Value Relation form controls.
- In its configuration, we find what is the last table (PESTICIDE), what field is displayed (NAME) and what field is used as ID.
The relation will be used inside ANALYSIS form to make the link between ANALYSIS and ANALYSIS_PESTICIDE. The QListWidget needs to have the name of the relation for the code to have a way to find which tables are involved. Value Relation is for the other part of the n,n relation: ANALYSIS_PESTICIDE to PESTICIDE.
We are done with the concepts !
Show me the code !
n,n dedicated dialog
Time to dive into Python...
First thing to do: the pesticide UI ! With PyQt you can create a .ui file and build it with qt4-designer. But loading a .ui file from Python can be unsafe: you have to deal with the file location. As the pesticide UI is very trivial, I prefer to build it with code. So, here is the Python code of the dialog:
def setupUi(self):
'''Builds the QDialog'''
# Form building
self.setObjectName(u"nnDialog")
self.resize(550, 535)
self.setMinimumSize(QtCore.QSize(0, 0))
self.buttonBox = QtGui.QDialogButtonBox(self)
self.buttonBox.setGeometry(QtCore.QRect(190, 500, 341, 32))
self.buttonBox.setOrientation(QtCore.Qt.Horizontal)
self.buttonBox.setStandardButtons(QtGui.QDialogButtonBox.Cancel|QtGui.QDialogButtonBox.Ok)
self.buttonBox.setObjectName(u"buttonBox")
self.verticalLayoutWidget = QtGui.QWidget(self)
self.verticalLayoutWidget.setGeometry(QtCore.QRect(9, 9, 521, 491))
self.verticalLayoutWidget.setObjectName(u"verticalLayoutWidget")
self.verticalLayout = QtGui.QVBoxLayout(self.verticalLayoutWidget)
self.verticalLayout.setMargin(0)
self.verticalLayout.setObjectName(u"verticalLayout")
self.horizontalLayout = QtGui.QHBoxLayout()
self.horizontalLayout.setObjectName(u"horizontalLayout")
self.label = QtGui.QLabel(self.verticalLayoutWidget)
self.label.setObjectName(u"label")
self.horizontalLayout.addWidget(self.label)
self.SEARCH = QtGui.QLineEdit(self.verticalLayoutWidget)
self.SEARCH.setObjectName(u"SEARCH")
self.horizontalLayout.addWidget(self.SEARCH)
self.verticalLayout.addLayout(self.horizontalLayout)
self.horizontalLayout_2 = QtGui.QHBoxLayout()
self.horizontalLayout_2.setObjectName(u"horizontalLayout_2")
self.LIST = QtGui.QListWidget(self.verticalLayoutWidget)
self.LIST.setObjectName(u"LIST")
self.horizontalLayout_2.addWidget(self.LIST)
self.verticalLayout.addLayout(self.horizontalLayout_2)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.rejected.connect(self.reject)
QtCore.QMetaObject.connectSlotsByName(self)
Well, this is what you can have from pyuic4 from a qt4-designer .ui file. But this time you don't need the file anymore.
n,n list behaviour functions
Next thing to do is to populate the n,n dialog. We also need to add the "search engine" functions and a way to pre-check values that are in the ANALYSIS_PESTICIDE table for the current analysis. And at the end, we need to send the checked values to the main form (ANALYSIS one) in order to make the database update and to update the form control.
I've created a class for this:
class nnDialog(QtGui.QDialog):
'''Dedicated n,n relations Form Class'''
def __init__(self, parent, layer, shownField, IdField, initValues, search=False):
'''Constructor'''
QtGui.QDialog.__init__(self,parent)
self.initValues = initValues
self.shownField = shownField
self.layer = layer
self.IdField = IdField
self.search = search
if self.layer is None and DEBUGMODE:
QgsMessageLog.logMessage(u"nnDialog constructor: The layer {0} doesn't exists !".format(layer.name()),"nnForms", QgsMessageLog.INFO)
# Build the GUI and populate the list with the good values
self.setupUi()
self.populateList()
# Add dynamic control when list is changing
self.SEARCH.textChanged.connect(self.populateList)
self.LIST.itemChanged.connect(self.changeValues)
def setupUi(self):
'''Builds the QDialog'''
# Form building
self.setObjectName(u"nnDialog")
self.resize(550, 535)
self.setMinimumSize(QtCore.QSize(0, 0))
self.buttonBox = QtGui.QDialogButtonBox(self)
self.buttonBox.setGeometry(QtCore.QRect(190, 500, 341, 32))
self.buttonBox.setOrientation(QtCore.Qt.Horizontal)
self.buttonBox.setStandardButtons(QtGui.QDialogButtonBox.Cancel|QtGui.QDialogButtonBox.Ok)
self.buttonBox.setObjectName(u"buttonBox")
self.verticalLayoutWidget = QtGui.QWidget(self)
self.verticalLayoutWidget.setGeometry(QtCore.QRect(9, 9, 521, 491))
self.verticalLayoutWidget.setObjectName(u"verticalLayoutWidget")
self.verticalLayout = QtGui.QVBoxLayout(self.verticalLayoutWidget)
self.verticalLayout.setMargin(0)
self.verticalLayout.setObjectName(u"verticalLayout")
self.horizontalLayout = QtGui.QHBoxLayout()
self.horizontalLayout.setObjectName(u"horizontalLayout")
self.label = QtGui.QLabel(self.verticalLayoutWidget)
self.label.setObjectName(u"label")
self.horizontalLayout.addWidget(self.label)
self.SEARCH = QtGui.QLineEdit(self.verticalLayoutWidget)
self.SEARCH.setObjectName(u"SEARCH")
self.horizontalLayout.addWidget(self.SEARCH)
self.verticalLayout.addLayout(self.horizontalLayout)
self.horizontalLayout_2 = QtGui.QHBoxLayout()
self.horizontalLayout_2.setObjectName(u"horizontalLayout_2")
self.LIST = QtGui.QListWidget(self.verticalLayoutWidget)
self.LIST.setObjectName(u"LIST")
self.horizontalLayout_2.addWidget(self.LIST)
self.verticalLayout.addLayout(self.horizontalLayout_2)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.rejected.connect(self.reject)
QtCore.QMetaObject.connectSlotsByName(self)
def changeValues(self, element):
'''Whenever a checkbox is checked, modify the values'''
# Check if we check or uncheck the value:
if element.checkState() == Qt.Checked:
self.initValues.append(element.data(Qt.UserRole))
else:
self.initValues.remove(element.data(Qt.UserRole))
def populateList(self, txtFilter=None):
'''Fill the QListWidget with values'''
# Delete everything
self.LIST.clear()
# We need a request
request = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry)
if txtFilter is not None:
fields = self.layer.dataProvider().fields()
fieldname = fields[self.shownField].name()
request.setFilterExpression(u"\"{0}\" LIKE '%{1}%'".format(fieldname, txtFilter))
# Grab the results from the layer
features = self.layer.getFeatures(request)
for feature in sorted(features, key = lambda f: f[0]):
attr = feature.attributes()
value = attr[self.shownField]
element = QListWidgetItem(value)
element.setData(Qt.UserRole, attr[self.IdField])
# initValues will be checked
if attr[self.IdField] in self.initValues:
element.setCheckState(Qt.Checked)
else:
element.setCheckState(Qt.Unchecked)
self.LIST.addItem(element)
def getValues(self):
'''Return the selected values of the QListWidget'''
return self.initValues
The class is named nnDialog and it deals with the n,n dialog used to add/remove pesticides of the current analysis. The constructor is very simple:
- We need to know which layer will be displayed,
- what is the name of the field that will be displayed in the list,
- what is the name of the field used as ID,
- what are the values already checked (stored into ANALYSIS_PESTICIDE)
- once everything is transmitted by arguments to the constructor, we have to create the UI (see above),
- populate the list
- and add dynamic controls for search QLineEdit and QListWidget.
The changeValues method is called when you check a checkBox in the list. Whenever there is action, the ID_PESTICIDE value is added/removed from initValues.
The populateList method is used when the n,n dialog is opened (called by the constructor) and whenever there is some changes in the search text bar. This method is used to populate the list:
- List is first cleared.
- If there is some text in the search QLineEdit, we make a request on the displayed field (NAME is our case) of the current layer (PESTICIDE) to retrieve only the correct values.
- Otherwise, we grab all the values of the layer.
- We sort them alphabetically.
- And for each value, we create a QListWidgetItem (element of a list) with a checkBox.
- If the value is in the initValues, it is checked, otherwise, it is unchecked.
nnDialog class implements all the logic of the n,n Dialog and it's code is quite generic: every parameters are transmitted by the constructor. This class is used when you click on the "Modify" button at the right of the list of pesticides in the ANALYSIS form.
But to stay generic we have also to be generic with the code which triggers nnDialog...
main form code
Last thing to do: add logic to the ANALYSIS form. Here is the code:
class nnForm:
'''Class to handle forms to type data'''
def __init__(self, dialog, layerid, featureid):
self.dialog = dialog
self.layerid = layerid
self.featureid = featureid
self.nullValue = QSettings().value("qgis/nullValue" , u"NULL")
self.search = False
def id2listWidget(self, table, values, listWidget):
'''Show all the selected values of a link table on a QListWidget'''
# Find the Widget
if listWidget is None or table is None:
QgsMessageLog.logMessage(u"id2listWidget: We need to have a relation and a true widget !", "DBPAT", QgsMessageLog.INFO)
return False
# Empty the list
listWidget.clear()
# Get the params (for the first child table)
if self.valueRelationParams(table):
params = self.valueRelationParams(table)[0]
if params is None or not params:
QgsMessageLog.logMessage(u"id2listWidget: You need to add Value Relation to layer: {0} !".format(table.name()), "nnForms", QgsMessageLog.INFO)
return False
# Get target layer:
tgtLayer = params['tgtLayer']
# Handle values: need to escape \' characters
values = [v.replace(u"'", u"''") if isinstance(v, basestring) else v for v in values]
## Then, get the real values from other-side table
if values:
request = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry)
if params[u'tgtIdType'] in (QVariant.String, QVariant.Char):
query = u"{0} IN ('{1}')".format(params[u'tgtId'], u"','".join(values))
else:
query = u"{0} IN ({1})".format(params[u'tgtId'], u",".join([unicode(x) for x in values]))
request.setFilterExpression(query)
# and display them in the QListWidget
for feature in tgtLayer.getFeatures(request):
value = feature.attributes()[params[u'tgtValueIdx']]
if value != u"NULL":
element = QListWidgetItem(value)
element.setData(Qt.UserRole, feature.attributes()[params[u'tgtIdIdx']])
listWidget.addItem(element)
return True
def valueRelationParams(self,layer):
'''Function that returns the configuration parameters of a valueRelation as a list of dict'''
params = []
if layer is not None:
for idx, field in enumerate(layer.dataProvider().fields()):
if layer.editorWidgetV2(idx) == u"ValueRelation":
param = {}
param[u'srcId'] = field.name()
param[u'srcIdIdx'] = idx
if u"Layer" in layer.editorWidgetV2Config(idx):
tgtLayerName = layer.editorWidgetV2Config(idx)[u"Layer"]
tgtLayer = QgsMapLayerRegistry.instance().mapLayer(tgtLayerName)
if tgtLayer is None:
QgsMessageLog.logMessage(u"valueRelationParams: Can't find the layer {0} !".format(tgtLayerName), "nnForms", QgsMessageLog.INFO)
return False
param[u'tgtLayer'] = tgtLayer
param[u'tgtId'] = layer.editorWidgetV2Config(idx)[u"Key"]
param[u'tgtValue'] = layer.editorWidgetV2Config(idx)[u"Value"]
# Find index of all fields:
for indx, f in enumerate(tgtLayer.dataProvider().fields()):
if f.name() == param[u'tgtId']:
param[u'tgtIdIdx'] = indx
param[u'tgtIdType'] = f.type()
if f.name() == param[u'tgtValue']:
param[u'tgtValueIdx'] = indx
params.append(param)
# notification
if not params:
QgsMessageLog.logMessage(u"valueRelationParams: There is not Value Relation for the layer {0} !".format(layer.name()), "nnForms", QgsMessageLog.INFO)
return params
def manageMultiple(self):
'''Handle specifics thesaurus form'''
# Scan all of the QgsRelations of the project
relations = QgsProject.instance().relationManager().relations()
for listWidget in [f for f in self.dialog.findChildren(QListWidget) if u"REL_" in f.objectName()]:
listName = listWidget.objectName()
if listName not in relations.keys():
QgsMessageLog.logMessage(u"manageMultiple: There is no Relation for control {0} !".format(listWidget.objectName()), "nnforms", QgsMessageLog.INFO)
continue
# Find what is the table to show
relation = relations[listName]
shownLayer = relation.referencingLayer()
# Find other side of n,n relation
if self.valueRelationParams(shownLayer):
params = self.valueRelationParams(shownLayer)[0]
if params is None:
continue
# When found, we are ready to populate the QListWidget with the good values
values = []
if self.featureid:
# Get the features to display
request = relation.getRelatedFeaturesRequest(self.featureid)
request.setFlags(QgsFeatureRequest.NoGeometry)
for feature in shownLayer.getFeatures(request):
values.append(feature.attributes()[params[u'srcIdIdx']])
self.id2listWidget(shownLayer, values, listWidget)
buttonWidget = self.dialog.findChild(QPushButton, listName+u"_B")
if buttonWidget:
if self.search or self.layerid.isEditable():
buttonWidget.clicked.connect(partial(self.openSubform, listWidget, relation, values))
buttonWidget.setEnabled(True)
else:
buttonWidget.setEnabled(False)
elif DEBUGMODE:
QgsMessageLog.logMessage(u"manageMultiple: There is no button for control {0} !".format(listName), "nnForms", QgsMessageLog.INFO)
def openSubform(self, widget, relation, values):
'''Open a dedicated dialog form with values taken from a child table.'''
table = relation.referencingLayer()
if self.valueRelationParams(table):
params = self.valueRelationParams(table)[0]
if params is None or not params:
QgsMessageLog.logMessage(u"openSubform: There is no Value Relation for layer: {0} !".format(table.name()), "nnForms", QgsMessageLog.INFO)
return False
if widget is None:
QgsMessageLog.logMessage(u"openSubForm: no widgets found for field {0} !".format(field), "nnForms", QgsMessageLog.INFO)
# Open the form with the good values
dialog = nnDialog(self.dialog, params[u'tgtLayer'], params[u'tgtValueIdx'], params[u'tgtIdIdx'], values, self.search)
# handle results
if dialog.exec_():
# Get the results:
thevalues = dialog.getValues()
# Modify target table if we have a featureid
if self.featureid:
table.startEditing()
caps = table.dataProvider().capabilities()
## Delete all the previous values
if caps & QgsVectorDataProvider.DeleteFeatures:
request = relation.getRelatedFeaturesRequest(self.featureid)
request.setFlags(QgsFeatureRequest.NoGeometry)
fids = [f.id() for f in table.getFeatures(request)]
table.dataProvider().deleteFeatures(fids)
## Add the new values
if caps & QgsVectorDataProvider.AddFeatures:
for value in thevalues:
feat = QgsFeature()
feat.setAttributes([None, self.featureid.attributes()[0], value])
table.dataProvider().addFeatures([feat])
## Commit changes
table.commitChanges()
# refresh listWidget aspect
self.id2listWidget(table, thevalues, widget)
The nnForm class will manage the form of ANALYSIS (or every form that has the same class Python function).
The manageMultiple method, will "scan" the layer form to find all QListWidgets with the same name than a relation. For each of those QListWidgets, we try to find what is the intermediate table (ANALYSIS_PESTICIDE) and what is the last table (from Value Relation). Then the QListWidget is populated with the values from ANALYSIS_PESTICIDE (and by retreiving the pesticides names). At last, the QPushButton that is named like the relation (+_B) is connected to a method which will open a nnDialog (see previous chapter).
OpenSubForm method is used to create the nnDialog (from the same named class), to give it the already checked values and to grab the result once the nnDialog dialog is closed. Most of the code of this method is for updating values with quite a brutal approach: we erase every data stored into ANALYSIS_PESTICIDE that have the same ID_ANALYSIS value than the current analysis ! Then, we re-add everything... But it seems to be faster than filtering the already checked values ! At last, th QListWidget involved is refreshed.
id2listWidget is the method used to populate and refresh a QListWidget with relations on the form. Everything is first cleared. A request to the last table is done (PESTICIDE) to grab the field that msut be shown (NAME). The values (IDs) are requested before and put into the constructor of this method.
valueRelationParams is used to find what are: the target layer, the shown field, the identifying field of a value relation control configuration of a table. It is used in manageMultiple and id2listWidget methods to find what to display.
Putting everything into one file
# -*- coding: utf-8 -*-
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import QgsMapLayerRegistry, QgsMessageLog, QgsFeatureRequest, QgsFeature
from qgis.core import QgsRelationManager, QgsRelation, QgsProject, QgsVectorDataProvider
from qgis.utils import iface
from functools import partial
from PyQt4 import QtCore, QtGui
# Global variables
DEBUGMODE = True
class nnDialog(QtGui.QDialog):
'''Dedicated n,n relations Form Class'''
def __init__(self, parent, layer, shownField, IdField, initValues, search=False):
'''Constructor'''
QtGui.QDialog.__init__(self,parent)
self.initValues = initValues
self.shownField = shownField
self.layer = layer
self.IdField = IdField
self.search = search
if self.layer is None and DEBUGMODE:
QgsMessageLog.logMessage(u"nnDialog constructor: The layer {0} doesn't exists !".format(layer.name()),"Your App", QgsMessageLog.INFO)
# Build the GUI and populate the list with the good values
self.setupUi()
self.populateList()
# Add dynamic control when list is changing
self.SEARCH.textChanged.connect(self.populateList)
self.LIST.itemChanged.connect(self.changeValues)
def setupUi(self):
'''Builds the QDialog'''
# Form building
self.setObjectName(u"nnDialog")
self.resize(550, 535)
self.setMinimumSize(QtCore.QSize(0, 0))
self.buttonBox = QtGui.QDialogButtonBox(self)
self.buttonBox.setGeometry(QtCore.QRect(190, 500, 341, 32))
self.buttonBox.setOrientation(QtCore.Qt.Horizontal)
self.buttonBox.setStandardButtons(QtGui.QDialogButtonBox.Cancel|QtGui.QDialogButtonBox.Ok)
self.buttonBox.setObjectName(u"buttonBox")
self.verticalLayoutWidget = QtGui.QWidget(self)
self.verticalLayoutWidget.setGeometry(QtCore.QRect(9, 9, 521, 491))
self.verticalLayoutWidget.setObjectName(u"verticalLayoutWidget")
self.verticalLayout = QtGui.QVBoxLayout(self.verticalLayoutWidget)
self.verticalLayout.setMargin(0)
self.verticalLayout.setObjectName(u"verticalLayout")
self.horizontalLayout = QtGui.QHBoxLayout()
self.horizontalLayout.setObjectName(u"horizontalLayout")
self.label = QtGui.QLabel(self.verticalLayoutWidget)
self.label.setObjectName(u"label")
self.horizontalLayout.addWidget(self.label)
self.SEARCH = QtGui.QLineEdit(self.verticalLayoutWidget)
self.SEARCH.setObjectName(u"SEARCH")
self.horizontalLayout.addWidget(self.SEARCH)
self.verticalLayout.addLayout(self.horizontalLayout)
self.horizontalLayout_2 = QtGui.QHBoxLayout()
self.horizontalLayout_2.setObjectName(u"horizontalLayout_2")
self.LIST = QtGui.QListWidget(self.verticalLayoutWidget)
self.LIST.setObjectName(u"LIST")
self.horizontalLayout_2.addWidget(self.LIST)
self.verticalLayout.addLayout(self.horizontalLayout_2)
self.buttonBox.accepted.connect(self.accept)
self.buttonBox.rejected.connect(self.reject)
QtCore.QMetaObject.connectSlotsByName(self)
def changeValues(self, element):
'''Whenever a checkbox is checked, modify the values'''
# Check if we check or uncheck the value:
if element.checkState() == Qt.Checked:
self.initValues.append(element.data(Qt.UserRole))
else:
self.initValues.remove(element.data(Qt.UserRole))
def populateList(self, txtFilter=None):
'''Fill the QListWidget with values'''
# Delete everything
self.LIST.clear()
# We need a request
request = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry)
if txtFilter is not None:
fields = self.layer.dataProvider().fields()
fieldname = fields[self.shownField].name()
request.setFilterExpression(u"\"{0}\" LIKE '%{1}%'".format(fieldname, txtFilter))
# Grab the results from the layer
features = self.layer.getFeatures(request)
for feature in sorted(features, key = lambda f: f[0]):
attr = feature.attributes()
value = attr[self.shownField]
element = QListWidgetItem(value)
element.setData(Qt.UserRole, attr[self.IdField])
# initValues will be checked
if attr[self.IdField] in self.initValues:
element.setCheckState(Qt.Checked)
else:
element.setCheckState(Qt.Unchecked)
self.LIST.addItem(element)
def getValues(self):
'''Return the selected values of the QListWidget'''
return self.initValues
class nnForm:
'''Class to handle forms to type data'''
def __init__(self, dialog, layerid, featureid):
self.dialog = dialog
self.layerid = layerid
self.featureid = featureid
self.nullValue = QSettings().value("qgis/nullValue" , u"NULL")
self.search = False
def id2listWidget(self, table, values, listWidget):
'''Show all the selected values of a link table on a QListWidget'''
# Find the Widget
if listWidget is None or table is None:
QgsMessageLog.logMessage(u"id2listWidget: We need to have a relation and a true widget !", "nnForms", QgsMessageLog.INFO)
return False
# Empty the list
listWidget.clear()
# Get the params (for the first child table)
if self.valueRelationParams(table):
params = self.valueRelationParams(table)[0]
if params is None or not params:
QgsMessageLog.logMessage(u"id2listWidget: You need to add Value Relation to layer: {0} !".format(table.name()), "nnForms", QgsMessageLog.INFO)
return False
# Get target layer:
tgtLayer = params['tgtLayer']
# Handle values: need to escape \' characters
values = [v.replace(u"'", u"''") if isinstance(v, basestring) else v for v in values]
## Then, get the real values from other-side table
if values:
request = QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry)
if params[u'tgtIdType'] in (QVariant.String, QVariant.Char):
query = u"{0} IN ('{1}')".format(params[u'tgtId'], u"','".join(values))
else:
query = u"{0} IN ({1})".format(params[u'tgtId'], u",".join([unicode(x) for x in values]))
request.setFilterExpression(query)
# and display them in the QListWidget
for feature in tgtLayer.getFeatures(request):
value = feature.attributes()[params[u'tgtValueIdx']]
if value != u"NULL":
element = QListWidgetItem(value)
element.setData(Qt.UserRole, feature.attributes()[params[u'tgtIdIdx']])
listWidget.addItem(element)
return True
def valueRelationParams(self,layer):
'''Function that returns the configuration parameters of a valueRelation as a list of dict'''
params = []
if layer is not None:
for idx, field in enumerate(layer.dataProvider().fields()):
if layer.editorWidgetV2(idx) == u"ValueRelation":
param = {}
param[u'srcId'] = field.name()
param[u'srcIdIdx'] = idx
if u"Layer" in layer.editorWidgetV2Config(idx):
tgtLayerName = layer.editorWidgetV2Config(idx)[u"Layer"]
tgtLayer = QgsMapLayerRegistry.instance().mapLayer(tgtLayerName)
if tgtLayer is None:
QgsMessageLog.logMessage(u"valueRelationParams: Can't find the layer {0} !".format(tgtLayerName), "nnForms", QgsMessageLog.INFO)
return False
param[u'tgtLayer'] = tgtLayer
param[u'tgtId'] = layer.editorWidgetV2Config(idx)[u"Key"]
param[u'tgtValue'] = layer.editorWidgetV2Config(idx)[u"Value"]
# Find index of all fields:
for indx, f in enumerate(tgtLayer.dataProvider().fields()):
if f.name() == param[u'tgtId']:
param[u'tgtIdIdx'] = indx
param[u'tgtIdType'] = f.type()
if f.name() == param[u'tgtValue']:
param[u'tgtValueIdx'] = indx
params.append(param)
# notification
if not params:
QgsMessageLog.logMessage(u"valueRelationParams: There is not Value Relation for the layer {0} !".format(layer.name()), "nnForms", QgsMessageLog.INFO)
return params
def manageMultiple(self):
'''Handle specifics thesaurus form'''
# Scan all of the QgsRelations of the project
relations = QgsProject.instance().relationManager().relations()
for listWidget in [f for f in self.dialog.findChildren(QListWidget) if u"REL_" in f.objectName()]:
listName = listWidget.objectName()
if listName not in relations.keys():
QgsMessageLog.logMessage(u"manageMultiple: There is no Relation for control {0} !".format(listWidget.objectName()), "nnforms", QgsMessageLog.INFO)
continue
# Find what is the table to show
relation = relations[listName]
shownLayer = relation.referencingLayer()
# Find other side of n,n relation
if self.valueRelationParams(shownLayer):
params = self.valueRelationParams(shownLayer)[0]
if params is None:
continue
# When found, we are ready to populate the QListWidget with the good values
values = []
if self.featureid:
# Get the features to display
request = relation.getRelatedFeaturesRequest(self.featureid)
request.setFlags(QgsFeatureRequest.NoGeometry)
for feature in shownLayer.getFeatures(request):
values.append(feature.attributes()[params[u'srcIdIdx']])
self.id2listWidget(shownLayer, values, listWidget)
buttonWidget = self.dialog.findChild(QPushButton, listName+u"_B")
if buttonWidget:
if self.search or self.layerid.isEditable():
buttonWidget.clicked.connect(partial(self.openSubform, listWidget, relation, values))
buttonWidget.setEnabled(True)
else:
buttonWidget.setEnabled(False)
elif DEBUGMODE:
QgsMessageLog.logMessage(u"manageMultiple: There is no button for control {0} !".format(listName), "nnForms", QgsMessageLog.INFO)
def openSubform(self, widget, relation, values):
'''Open a dedicated dialog form with values taken from a child table.'''
table = relation.referencingLayer()
if self.valueRelationParams(table):
params = self.valueRelationParams(table)[0]
if params is None or not params:
QgsMessageLog.logMessage(u"openSubform: There is no Value Relation for layer: {0} !".format(table.name()), "nnForms", QgsMessageLog.INFO)
return False
if widget is None:
QgsMessageLog.logMessage(u"openSubForm: no widgets found for field {0} !".format(field), "nnForms", QgsMessageLog.INFO)
# Open the form with the good values
dialog = nnDialog(self.dialog, params[u'tgtLayer'], params[u'tgtValueIdx'], params[u'tgtIdIdx'], values, self.search)
# handle results
if dialog.exec_():
# Get the results:
thevalues = dialog.getValues()
# Modify target table if we have a featureid
if self.featureid:
table.startEditing()
caps = table.dataProvider().capabilities()
## Delete all the previous values
if caps & QgsVectorDataProvider.DeleteFeatures:
request = relation.getRelatedFeaturesRequest(self.featureid)
request.setFlags(QgsFeatureRequest.NoGeometry)
fids = [f.id() for f in table.getFeatures(request)]
table.dataProvider().deleteFeatures(fids)
## Add the new values
if caps & QgsVectorDataProvider.AddFeatures:
for value in thevalues:
feat = QgsFeature()
feat.setAttributes([None, self.featureid.attributes()[0], value])
table.dataProvider().addFeatures([feat])
## Commit changes
table.commitChanges()
# refresh listWidget aspect
self.id2listWidget(table, thevalues, widget)
def opennnForm(dialog, layerid, featureid):
'''Generic function to open a nnForm'''
form = nnForm(dialog, layerid, featureid)
QgsMessageLog.logMessage(u"opennnForm !", "nnforms", QgsMessageLog.INFO)
form.manageMultiple()
Conclusion
Okay, this one is quite complex ! If you want to implement n,n forms, you have to code because QGIS is not able to handle them for the moment. For a true implementation into QGIS code, I would take a different path. I can imagine to have a new relation type dedicated to n,n relations. In those relations, you would have to:
- Declare the "parent" layer (ANALYSIS in our case).
- Declare the "intermediate" table (ANALYSIS_PESTICIDE) and the shared attributes.
- Declare the "displayed" layer (PESTICIDE) and the shared attributes (with the intermediate table).
Once in the form, you would use a new form control to configure:
- the displayed field(s) or expression(s).
- if you want a search bar or not.
- if you want to filter values of the "displayed" layer.
Introduction
This article is the next of my serie of articles about QGIS forms. You can read the previous one here. Today we will focus on photo form controls. In QGis, they are dedicated to image viewing. Change your Edit Widget to Photo. You can then edit the size of the photo display (seems to be required if you want to effectively display the photo on the form).
Nice widget, isn't it ? Be we can go further...
Custom UI
With this configuration, QGis is able to render photos on your Widget by its own means. You do not need to code. But you perhaps want to go further and show some url links for documents that are not photos or images but other things (think about a photo inside a PDF file or a video) ? Or perhaps, you would like to use a Photo widget into a custom UI ? Let's see how it can be achieved.
The first thing to do is to build a custom form with qtcreator (designer-qt4 i Debian Jessie). You'll have to create a QWidget. It will be named with the field you have used to store the path to the photo in the data layer. This QWidget will hold a QGridLayout. Inside this layout, put a QLabel and name it PhotoLabel (required name). It will be used to display the photo. Then, add a QLineEdit named lineEdit (required name). It will be used to show the path of the photo. You also need to add a QPushButton named FileChooserButton in the layout.
If you want to show an URL, add a QLabel and name it with the name of your choice. You'll have to be careful to add two properties for this QLabel:
- openExternalLinks: checked (otherwise, you will not be able to open the file with a click).
- textInteractionFlags: LinksAccessibleByMouse.
With this form, QGIS will be able to display the photo in the fields and the widget will work as if it was auto-generated by QGIS.
Code
Time to dive into Python... We want to show a URL to open an external document. This document will be our image (that is presently shown) or another type of document (PDF/video/sound/etc.).
def manageIll(dialog, layerid, featureid):
'''Handle link for files in Illustration'''
# Find file value:
child = dialog.findChild(QLineEdit, u"lineEdit")
if child:
child.textChanged.connect(partial(modifyPhoto, dialog))
modifyPhoto(dialog)
def modifyPhoto(dialog):
'''Function to buidl link for files in Illustration form'''
# Find the fie path value
lineEdit = dialog.findChild(QLineEdit, u"lineEdit")
if lineEdit:
fileValue = lineEdit.text()
nullValue = QSettings().value("qgis/nullValue" , u"NULL")
if fileValue == nullValue or fileValue is None:
return False
basename = os.path.basename(fileValue)
filename = u"<a href=\"file:///{0}\">{1}</a>".format(fileValue,basename)
# Affect the value to the URL QLabel
urlLabel = dialog.findChild(QLabel, u"ILLURL_L")
if urlLabel and fileValue is not None:
urlLabel.setText(filename)
# Determine if we can produce a QPixmap from the file
if not QPixmap().load(fileValue):
photoLabel = dialog.findChild(QLabel, u"PhotoLabel")
if photoLabel:
photoLabel.setText(u"This file is not a photo !")
The manageIll function will connect the signal textChanged of the QLineEdit named "lineEdit" (remember the UI part) to a dedicated function named modifyPhoto. modifyPhoto will take the value of "lineEdit" and will transform it into an "URL" that will be shown into a QLabel named ILLURL_L (see UI part). At the end, we try to find if the filepath from "lineEdit" is a valid image and if it is not the case, we display a text inside the dedicated photo QLabel (named PhotoLabel).
You'll need to change the objectName property to reflect the widgets names found in your form (and in your layer of course).
Conclusion
This article just show an easy way to add a QGis photo form control inside a custom form. You can also add code to change the default behaviour of your photo control with Python. Whitout doubt, I am sure that one day, the QGIS developpers will implement a better external documents behaviour in QGIS...
Introduction
A la suite de mes différents articles sur l'installation et l'empaquetage du client Oracle, il est temps de passer à l'objectif final de ces opérations: installer le provider Oracle Spatial (que j'appelle aussi connecteur Oracle Spatial) pour QGIS sous Debian Jessie.
Nous allons récupérer les sources de QGIS, basculer sur une branche stable, effectuer nos petites modifications pour inclure les bibliothèques Oracle et créer les différents paquets. A l'issue de cette procédure, nous disposerons d'une installation complète de QGIS avec la possibilité de se connecter à un serveur Oracle Spatial, que ce soit pour un usage station de travail ou pour utiliser QGIS Server.
Mode opératoire
Pour ma part, j'effectue l'empaquetage dans une machine virtuelle. En effet, cette opération d'empaquetage de QGIS nécéssite d'installer une sacrée liste de paquets et je ne souhaite pas alourdir ma station de travail avec des paquets qui ne me serviront qu'une seule fois (ou dans tous les cas assez peu souvent).
Récupérer les sources et pointer vers la bonne version de QGIS
$ mkdir -p ~/packaging
$ git clone https://github.com/qgis/QGIS.git
$ cd ~/packaging
$ git clone ../QGIS
$ git checkout final-2_10_0
C'est la version 2.10 qui sera utilisée dans la suite de cet article mais les instructions doivent fonctionner aussi bien pour les versions ultérieures que pour les versions antérieures, à partir de la version 2.8.
Paquets à installer
Avant de pouvoir empaqueter QGIS, il va vous falloir un sacré paquet de paquets !
Voici la liste:
# apt install build-essential ca-certificates devscripts fakeroot bison cmake debhelper flex grass-dev libexpat1-dev libfcgi-dev libgdal-dev libgeos-dev libgsl0-dev libpq-dev libproj-dev libqt4-dev libqt4-opengl-dev libqtwebkit-dev libqwt-dev libspatialite-dev libsqlite3-dev libspatialindex-dev pkg-config pyqt4-dev-tools python-all-dev python-dev python-qt4-dev python-sip-dev txt2tags doxygen python-qscintilla2 pyqt4.qsci-dev libosgearth-dev libopenscenegraph-dev libqscintilla2-dev graphviz xvfb xauth xfonts-base xfonts-100dpi xfonts-75dpi xfonts-scalable spawn-fcgi lighttpd poppler-utils python-pyspatialite qt4-doc-html libqt4-sql-sqlite python-matplotlib osgearth-data python-psycopg2 python-httplib2 python-jinja2 python-markupsafe liblwgeom-2.1.7 qt4-designer
Si vous avez suivi mon précédent article, vous pourrez bien entendu, installer les paquets oracle-instantclient et oracle-instantclient-dev pour pouvoir compiler le provider Oracle. Vous pouvez également simplement installer le client Oracle à la main.
Modification du fichier debian/control.in
Le fichier debian/control.in est un fichier modèle pour le fichier debian/control. Pour ceux qui font régulièrement de la fabrication de paquet Debian, ce fichier ne devrait pas avoir de secret. QGIS présente la particularité de mettre à part la construction du paquet dédié au connecteur Oracle Spatial. En effet, le client et le SDK Oracle étant non libres, ils ne sont pas empaquetés dans la distribution Debian.
Vous devez donc modifier le fichier debian/control.in pour générer le paquet du provider Oracle. Si vous l'omettez, le paquet ne sera pas construit et le connecteur ne sera pas disponible même si l'iĉone d'accès à la boîte de dialogue de sélection des couches sera présente dans l'interface de QGIS. Ajoutez les éléments à la fin du fichier:
Package: qgis-oracle-provider
Architecture: any
Depends: ${shlibs:Depends}, ${misc:Depends}, oracle-instantclient
Section: contrib/database
Description: QGIS oracle provider
QGIS is a Geographic Information System (GIS) which manages, analyzes and
display databases of geographic information.
.
This package contains the QGIS oracle provider.
Pensez à supprimer la dépendance vers oracle-instantclient si vous n'avez pas installé le client Oracle via le paquet non officiel que j'ai présenté dans cet article.
Modification du fichier debian/rules
Une fois le paquet du connecteur Oracle correctement déclaré, il reste à modifier le fichier debian/rules. En effet, par défaut, le fichier rules indique de ne fabriquer le paquet que sous la distribution sid-oracle. Or, nous sommes sous Jessie. Par défaut, la construction de paquet ne saura pas fabriquer le connecteur Oracle.
Vous devez modifier la variale CMAKE_OPTS et y inclure la directive -DWITH_ORACLE=TRUE, comme dans ce qui suit:
...
CMAKE_OPTS := \
-DBUILDNAME=$(DEB_BUILD_NAME) \
-DCMAKE_VERBOSE_MAKEFILE=1 \
-DCMAKE_INSTALL_PREFIX=/usr \
-DGRASS_PREFIX=/usr/lib/$(GRASS) \
-DBINDINGS_GLOBAL_INSTALL=TRUE \
-DPEDANTIC=TRUE \
-DWITH_QSPATIALITE=TRUE \
-DWITH_SERVER=TRUE \
-DWITH_SERVER_PLUGINS=TRUE \
-DSERVER_SKIP_ECW=TRUE \
-DQGIS_CGIBIN_SUBDIR=/usr/lib/cgi-bin \
-DWITH_APIDOC=TRUE \
-DWITH_CUSTOM_WIDGETS=TRUE \
-DWITH_ORACLE=TRUE \
-DWITH_GLOBE=TRUE \
-DWITH_INTERNAL_HTTPLIB2=FALSE \
-DWITH_INTERNAL_JINJA2=FALSE \
-DWITH_INTERNAL_MARKUPSAFE=FALSE \
-DWITH_INTERNAL_PYGMENTS=FALSE \
-DWITH_INTERNAL_DATEUTIL=FALSE \
-DWITH_INTERNAL_PYTZ=FALSE \
-DWITH_INTERNAL_SIX=FALSE
...
Modification de la configuration de CMake
Il faut modifier le fichier src/providers/oracle/ocispatial/cmake/FindOCI.cmake pour indiquer le répertoire du client Oracle et de son SDK. Si vous utilisez les paquets Debian du client Oracle OCI dont j'ai décris la réalisation, vous n'avez rien à faire. Si vous avez installé manuellement le client Oracle dans /opt, voici comment vous devez modifier le fichier:
FIND_PATH(OCI_INCLUDE_DIR oci.h
PATHS
/opt/oracle/instantclient_11_2/sdk/include
$ENV{OSGEO4W_ROOT}/include
$ENV{ORACLE_HOME}/rdbms/public
)
FIND_LIBRARY(OCI_LIBRARY clntsh oci
PATHS
/opt/oracle/instantclient_11_2/
$ENV{OSGEO4W_ROOT}/lib
$ENV{ORACLE_HOME}/lib
)
...
Création des paquets
Vous avez fait le plus dur ! Il ne reste qu'à lancer dpkg-buildpackage -us -uc -b et attendre ! Je vous conseille plutôt de lancer la commande qui suit:
$ DEB_BUILD_OPTIONS=nocheck dpkg-buildpackage -us -uc -b
L'option nocheck permet de ne pas jouer la batterie de tests à la suite de la compilation. Sachez que les tests prennent facilement une bonne demi-heure... QGIS est assez long à empaqueter comme ça (compter au moins une demi-heure sur une machine récente) !
Si jamais le processus échoue à un moment donné, il vous reste à "hacker" dans le répertoire debian et ailleurs pour voir ce qui manque. En règle générale, il s'agit d'un problème de dépendances. Pour éviter que dpkg-buildpackage ne relance toute la compilation from scratch, utilisez l'option -nc qui permet de ne pas nettoyer l'arborescence des sources. c'est très utile si vous avez un problème au moment de la création du paquet avec dh_shlibdeps par exemple.
A la fin du processus et si tout se passe bien, vous aurez une liste impressionnante de paquets deb dans le répertoire ~/packaging dont le fameux qgis-oracle-provider_2.10.0_amd64.deb !
Pour installer en mode brutal:
dpkg -i ~/packaging/*.deb
Tests
Lancez le client QGIS et tentez de créer une connexion Oracle Spatial. Si tout se passe bien, vous aurez accès à la liste de vos couches et vous pourrez en ajouter une dans le canevas de cartes.
Si ce n'est pas le cas, il va vous falloir "hacker" un peu plus. Le principal coupable sera sans doute le fichier tnsnames.ora. Lors de mon précédent article, j'avais indiqué qu'il fallait paramétrer une variable d'environnement TNS_ADMIN. Mais cette dernière n'est disponible que via un appel à Bash. Si vous lancez QGIS depuis le bureau (via les suckless-tools bien sûr), cette variable n'est pas initialisée et aucune connexion ne fonctionnera.
Vérifiez votre fichier ~/.profile pour voir s'il contient bien la variable TNS_ADMIN:
export TNS_ADMIN=~/.oracle
Une dernière recommandation
Pour assurer de bonnes performances sous Oracle Spatial, utilisez absolument les index spatiaux. Contrairement à PostGIS, les index spatiaux sont pratiquement des pré-requis sous Oracle Spatial car un très grand nombre de fonctions ne sont pas utilisables sans.
En plus de la création de l'index, vous devez vous assurer que l'utilisateur qui accède à vos données disposent également des droits de lecture sur les tables d'index (tables dont le nom commence par MDRT) sinon le provider Oracle de QGIS ne pourra pas utiliser les fonctions d'index. Voici un petit script en PL/SQL pour vous assurer de ces droits une fois tout vos index créés (remplacez utilisateur par le nom du compte Oracle à qui vous voulez donner accès à vos index, ce script doit être lancé par l'utilisateur propriétaire du schéma Oracle):
set serveroutput on
DECLARE
tbName VARCHAR2(200);
CURSOR mdrtTables IS
SELECT SDO_INDEX_TABLE FROM ALL_SDO_INDEX_METADATA;
BEGIN
-- We need to scan all the Index tables
FOR tb in mdrtTables
LOOP
DBMS_OUTPUT.PUT_LINE('Table: ' || tb.SDO_INDEX_TABLE);
EXECUTE IMMEDIATE 'GRANT SELECT ON ' || tb.SDO_INDEX_TABLE || ' TO utilisateur';
END LOOP;
END;
/
EXIT
Conclusion
Grâce à l'empaquetage Debian, il est possible d'inclure le connecteur Oracle de QGIS en modifiant quelques fichiers. Au final, le processus n'est pas si simple car il faut installer le client Oracle et recréer le paquet Debian à la main alors que si vous aviez un serveur PostgreSQL/PostGIS, tout serait prêt "out-of-the-box".
Mais bon, parfois, on n'a pas le choix ! J'ai bien galéré à trouver une méthode suffisamment propre et automatisée pour ce travail étant donné la complexité du projet QGIS et surtout la difficulté amenée par le provider Oracle qui s'appuie sur du logiciel non libre, mal intégré sous Debian.
Néanmoins, les entités qui souhaiteraient installer QGIS Server sous Debian et qui n'ont que des données géographiques disponibles sous Oracle pourront aller un peu plus loin grâce à ces instructions...
Introduction
Today I am going to inaugurate a new serie of articles about QGIS forms. QGIS is perhaps the best GIS software that ease the most attributes edition. You can build nearly complete GIS applications just by triggering some parameters on the QGIS fields dialog (and with a little bit of code, of course). I've made a whole application with QGIS forms. On the GIS side, it was not so hard: 8 layers to edit with simple geometries. But on the attributes side, it was a real challenge. Some facts about the database and the application:
- About one hundred attributes tables stored on Oracle Database.
- Some layers can have about 80 form controls.
- n,n relations.
- Specific tables to store tree data (more on this later).
- Complete custom form controls.
- We needed subforms to handle photos.
- Of course you need a true search engine on all of the attribute fields.
Building such an application was not very easy even if you already know QGIS well. At the beginning of the project I was not really sure that QGIS could match the trick... ...But at last, I just want to say that QGis works perfectly well on this application even with a lot of data loaded.
When I started to implement features into this QGis project, I often faced QGIS lack of features on the forms. I've made some bug or features report on http://hub.qgis.org. In this serie of articles, I will try to show you how to circumvent the different problems I faced.
About QGIS forms
Before diving into code, let's have some basic informations about QGIS forms...
Forms are built by QGIS in order to edit (non-geometric) attributes for an object on a defined layer. You can read more about the process in the official QGis documentation.
You need to understand how forms are working in QGIS. To my mind, there are three parts on this subject:
- Data fields: the goal of the form is to edit the values of the attributes in the layer. Data fields are named with the attributes name of the layer.
- GUI: This is the part that QGIS show to the user to edit alphanumerical data. You can choose between three types: auto-generated (QGIS build everything, for basic forms), drag'n drop designed (you build the form using tabs and groupboxes inside QGIS) or customised (you need to use QtCreator to build a .ui file that will be shown by QGIS). To make link between form controls and fields values, there is a trivial rule: form controls are named like the field they represent.
- Python code: QGIS allows you to add Python logic to the GUI part of the form. Whenever the form opens, you can call a Python function. This function can access to GUI objects using Qt functions but you also have the full power of QGIS Python API.
For Data fields definition, you have nothing to do: QGIS will use the layer definition to name attributes. For GUI, you have to configure (per-layer) the method (auto-generated/drag'n drop/custom) on the Fields tab of the layer properties dialog. You can read this introduction article to have further details...
For Python code, you just have to name the Python file (module) that will be used and the function that will be launched. Read this reference article to dig a little bit more.
Show a clickable URL on the form
Introduction
In QGIS forms, you can specify that a field will store a path to a file. You just have to use "File Name" for the Edit Widget. There is no option for this type of widget:
When you create a new feature on the layer, you will see the following control:
Push the "..." button and you open a file choosing dialog (like every Qt file dialog). When you have chosen the file, the path (/home/medspx/clint_eastwood.ods) is stored in the line edit (a QLineEdit) right to the field name (DOCUMENT).
But what if you want to open the file and not just store the path ? You have to trick QGIS a little bit.
Build a custom form like this
First of all you have to use custom form. Otherwise, QGis is using a QLineEdit to show the file path. We need to use a QLabel as only QLabel are able to show a URL link. In your .ui file just add the following to your layout:
As you can see, we have a QLabel to name the field ("Document"). Then you have another QLabel which contains "No file selected". Name this QLabel with a dedicated name (the objectName property will be DOCUMENT_URL). This QLabel will make the job of link URL presentation.Be careful to have a field which is named like the objectName otherwise, QGIS will not save the data into the layer.
You have to check two properties for this QLabel:
- openExternalLinks: checked (otherwise, you will not be able to open the file with a click).
- textInteractionFlags: LinksAccessibleByMouse.
Then we have the QPushButton which will be named DOCUMENT_URL_B. And that's all for the ui part. Now, it is time to dive into Python...
Code to show the file path as clickable link
First, the code:
def manageURL(dialog, layerid=None, featureid=None):
'''General function to manage URL in subforms'''
# Manage URL buttons
for child in [f for f in dialog.findChildren(QPushButton) if u"URL_B" in f.objectName()]:
try:
child.clicked.disconnect()
except:
pass
child.clicked.connect(partial(chooseFile, dialog, child.objectName()[:-2]))
# Make URL fields always clickable (even in non-edit mode)
for child in [f for f in dialog.findChildren(QLabel) if f.objectName()[-3:] == u"URL"]:
child.setEnabled(True)
def chooseFile(dialog, field):
'''Open a dedicated file picker form'''
lineEdit = dialog.findChild(QLabel, field)
if lineEdit is None:
QgsMessageLog.logMessage(u"chooseFile: There is no QLabel for field {0} !".format(field), "DBPAT", QgsMessageLog.INFO)
filename = QtGui.QFileDialog.getOpenFileName(dialog, u'Choisir un fichier...')
if filename is not None:
# parse the result to make a link
basename = os.path.basename(filename)
link = u"<a href=\"file:///{0}\">{1}</a>".format(filename,basename)
lineEdit.setText(link)
When you want to display links, you just have to use manageURL function to do the job.
This function will try to find all of the QPushButton that are named like "SOMETHING_URL_B". Those buttons will be disconnected from their previous "clicked" signal. Then, we will re-connect the "clicked" signal to our "chooseFile" function. This function will do the job of building the link from the file path. manageURL ends with another trick: if you are not in edit mode, QLabel with links are not clickable. So we have to manually enable every QLabel that are links (ends with "URL" in our case).
chooseFile function opens a file chooser dialog (a standard QFileDialog from Qt). If the file path is valid, we extract the filename (the basename, only the last part of the path) from the file path. And then we build a link which is a simple <a>. At last, we update the QLabel with the link content.
Here is a rendered thing of the link:
You can click on the link and it will open the file with the default operating system defined viewer.
Conclusion
This first article just demonstrate that QGIS is so wide open to external code that you can do nearly everything you want with the forms. Embedding PyQt4 and giving access to the full QGIS API from Python is very interesting because you can go further than standard QGIS forms without recompiling QGIS.
As you can see, once you have written a dedicated function to handle URL links for File controls, it is very easy to expand it to any layer you want: just correctly name a form control on your custom form file and add the correct path to the function.
In the following article, we will focus on Photo edit widget and we will try to use it in a custom form with a little bit of code as an extra...
Introduction
Since QGis 1.8 (and even before), you can use QGis to build forms to deal with the attributes when you create a new geographic object. Over the successive version of QGis, form building get plent of new features like a 1,1 value relation (v2.2) and now 1,n relations (v2.2). You can also use the autobuild feature with which QGis builds the form depending on what you put in the fields properties. There is alos a way to build you own form control with QTCreator ui files. Of course, you can control the form logic with Python code.
In the development of a specific form I had to build a cascading control: the value of one control is linked to the value of an upper control (controls are linked together). This is especially convenient for narrowing a choice. For example you want to choose a district or an area that restricts a list of counties that itself restricts a list of towns. I call this cascading form controls.
For the moment, QGis don't have an internal mechanism to deal with those controls. I've sent a feature request to implement it but I also found a workaround by using Python to control the controls content. Here I publish this solution.
Data presentation
I've got the following tables:
- On geographic table named OBJECTGEO.
- Flat tables
- One for AREAs
- One for Couties
- Last one for towns
The flat tables are linked together: each county has got a reference to an area (with an ID) and each town has got a reference to a county.
Here the schemas of all of the tables
Build the form
You can build the form on OBJECTGEO table just by following what is written below:
Be sure to add value relations controls for the ID_AREA, ID_COUNTY and ID_TOWN fields:
For the moment, your form looks like this:
You can find that the values of the form are not linked together (if you are French you know that Finistère is not in Pays-de-la-loire and that SAUMUR is not in Finistère). Time to change this...
Write some Python code
Now that you have some controls in your form, you need to add a little bit of logic in order to have cascading controls. Everything is made in the following Python code.
# -*- encoding: utf-8 -*-
from PyQt4.QtCore import *
from PyQt4.QtGui import *
myDialog = None
def formOpen(dialog,layerid,featureid):
global myDialog
myDialog = dialog
# We need to introspect the form
# In QGis v2.4 autobuild form controls do not have a name... (fixed in v2.6)
# So we need to retrieve form controls by field order !
lst_children = dialog.findChildren(QComboBox)
area = lst_children[0]
county = lst_children[1]
town = lst_children[2]
# Clear the children QComboBox widgets
county.clear()
town.clear()
# When you change the value for area widget, change the content of the county widget
area.currentIndexChanged.connect(lambda: cascadeManage('COUNTY', 'ID_COUNTY', 'COUNTY', 'ID_AREA', county, area))
# When you change the value for county widget, change the content of the town widget
county.currentIndexChanged.connect(lambda: cascadeManage('TOWN', 'ID_TOWN', 'TOWN', 'ID_COUNTY', town, county))
def cascadeManage(layername, id_idx, txt_idx, fltr, widget, parent_widget):
'''Generic function to manage cascading QComboBox.
* layername: the layer (flat table) to request values for this control
* id_idx: primary key for the table (used in for the form control key)
* txt_idx: field shown in the form control QComboBox
* fltr: field which is used to filter the content of this widget based on the parent widget value
* widget: widget object (QComboBox) to control
* parent_widget: get the value from the parent QComboBox widget'''
from qgis.core import QgsMapLayerRegistry, QgsExpression
# Get the layer (flat table)
layers = QgsMapLayerRegistry.instance().mapLayersByName(layername)
if (len(layers) > 0):
layer = layers[0]
else:
return False
# get attributes indexes
ididx = layer.fieldNameIndex(id_idx)
txtidx = layer.fieldNameIndex(txt_idx)
# Get the value of the parent QComboBox
filter_code = str(parent_widget.itemData(parent_widget.currentIndex()))
# Build an expression with it
exp = QgsExpression("{0} = {1}".format(fltr, filter_code))
exp.prepare(layer.pendingFields())
# clear the current widget content
widget.clear()
# And fill it with filtered data:
for feature in layer.getFeatures():
value = exp.evaluate(feature)
if exp.hasEvalError():
raise ValueError(exp.evalErrorString())
if bool(value):
widget.addItem(feature.attributes()[txtidx], feature.attributes()[ididx])
You can find explanations directly in the code...
In action
Save the upper code in a file named MYFORM.py in the same directory than your QGis project.
Now, make the link between the form and the code: change the value of Python Init function text with MYFORM.formOpen:
- MYFORM refers to MYFORM.py Python file.
- formOpen is the name of the method which will be launched when the form opens.
Now, you can have actives cascading controls:
Whenever you choose an area, it restricts the list of counties. Same thing with counties and towns...
Beware, there is a bug in QGis 2.4 which just duplicate the new object when you use a form with Python code (you can find twice on the attribute table). But it is already fixed in the master version of QGis which will become v2.6 in november 2014. For QGis 2.4, you will just delete the duplicate entry in the attribute table.
Conclusion
As a conclusion, we have a mean to implement cascading form controls even (and especially) with value-relation controls. The code could be greatly improved. For example, with the 2.6 version of QGis, you can use field name instead of trying to find the good widgets sorted by order. You could go a bit deeper and make your own ui form (with QtCreator) and create cascading form controls that only modify the value of one field (in our case, we add AREA and COUNTY fields in OBJECTGEO but the real information we want is the ID of the TOWN).
I hope that QGis will provide soon an intern (and easier to use without Python code) method to use cascading form controls !
Introduction
Since QGis 2.0, you can use GeoAlgorithms (aka Processing or geoprocessing) to manipulate your data. But you can also developp your own GeoAlgo to do things that are not (already) included in QGis. As GeoAlgo is relatively young, the development on this part of QGis is important.
Today I propose a very basic GeoAlgo that takes an Oracle layer, use internal Oracle Spatial procedure (SDO_GEOM.VALIDATE_GEOMETRY_WITH_CONTEXT) to find non valid geometries and returns a layer with those geometries and an explanation of the errors that have been found.
Install cx_Oracle
The code needs cx_Oracle Python library to work. This library is a low-level communication Python (DBAPIv2) library for interconnections to Oracle databases. For the moment it is not included in QGis binary installation, so you have to install by yourself first.
If you are running a GNU/Linux distribution and have installed QGis from the official repositories, you'll have to install Oracle Client and SDK and then install cx_Oracle. Read this...
If you are running QGis on MS-Windows, you have to download the good version of cx_Oracle binaries (32 or 64 bits) and extract and copy cx_Oracle.pyd file in C:\Program Files\QGis Chugiak/apps/Python2.7/DLL .
To check for a valid installation of cx_Oracle in the QGis Python environment:
- Launch QGis
- Open the Python Console (Plugins -> Python Console)
- Type:
import cx_Oracle
- If you have no error message, you're done !
Source code for VerifyOracle.py
Here is the source of the GeoAlgo. Just save it in a file named VerifyOracle.py ans place it in the good directory. You'll have to find the directory that stores the official Processing plugin.
On MS-Windows systems, it is stored in C:\Program Files\QGIS Chugiak\apps\qgis\python\plugins\processing\algs\qgis.
# -*- coding: utf-8 -*-
#***************************************************************************
# VerifyOracle.py
# ---------------------
# Date : September 2014
# Copyright : (C) 2014 by Médéric RIBREUX
# Email : mederic.ribreux@medspx.fr
#***************************************************************************
#* *
#* This program is free software; you can redistribute it and/or modify *
#* it under the terms of the GNU General Public License as published by *
#* the Free Software Foundation; either version 2 of the License, or *
#* (at your option) any later version. *
#* *
#***************************************************************************
#
#You need to install cx_Oracle under QGis Python directory.
__author__ = 'Médéric RIBREUX'
__date__ = 'September 2014'
__copyright__ = '(C) 2014, Médéric RIBREUX'
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'
from osgeo import gdal
import cx_Oracle
import re
from qgis.core import *
from PyQt4.QtCore import *
from processing.core.GeoAlgorithm import GeoAlgorithm
from processing.parameters.ParameterVector import ParameterVector
from processing.parameters.ParameterTableField import ParameterTableField
from processing.outputs.OutputVector import OutputVector
from processing.tools import dataobjects, vector
from processing.tools.general import *
class VerifyOracle(GeoAlgorithm):
INPUT_VECTOR = 'INPUT_VECTOR'
OUTPUT = 'OUTPUT'
FIELD = 'FIELD'
def defineCharacteristics(self):
self.name = 'Verify layer geometries with Oracle Spatial engine'
self.group = 'Vector analysis tools'
self.addParameter(ParameterVector(self.INPUT_VECTOR, 'Vector layer',
[ParameterVector.VECTOR_TYPE_ANY]))
self.addParameter(ParameterTableField(self.FIELD, 'UID Field for input vector layer',
self.INPUT_VECTOR))
self.addOutput(OutputVector(self.OUTPUT, 'Result Vector layer'))
def processAlgorithm(self, progress):
uri = self.getParameterValue(self.INPUT_VECTOR)
layer = dataobjects.getObjectFromUri(uri)
# Deals with fields
fieldName = self.getParameterValue(self.FIELD)
fieldIdx = layer.fieldNameIndex(fieldName)
fields = layer.dataProvider().fields()
# Add the Errors field
fields.append(QgsField('Errors', QVariant.String))
# Get connection parameters
regexp = re.compile(".*dbname='([^']+)'.*user='([^']+).*password='([^']+)'.*table=([^ ]+).*\(([^\)]+)\)")
dbname, user, password, table, geocol = regexp.match(uri).groups()
query = u"SELECT c.{0}, SDO_GEOM.VALIDATE_GEOMETRY_WITH_CONTEXT (c.{1}, 0.001) FROM {2} c WHERE SDO_GEOM.VALIDATE_GEOMETRY_WITH_CONTEXT (c.{1}, 0.001) <> 'TRUE'".format(self.getParameterValue(self.FIELD),geocol, table)
# Make the connection and the query
connection = cx_Oracle.connect( user, password, dbname)
c = connection.cursor()
c.execute(query)
rows = c.fetchall()
c.close()
# Open a writer to the output vector layer
writer = self.getOutputFromName(
self.OUTPUT).getVectorWriter(fields,
layer.dataProvider().geometryType(),
layer.dataProvider().crs())
# We get all of the features from the input vector layer
features = vector.features(layer)
# And make some computations for the progress bar
total = 100.0 / float(len(features))
current = 0
outFeat = QgsFeature()
# Build a list of errors (at least as big as number of features of the layer)
errors = []
for row in rows:
errors.append({'GID':row[0], 'ERROR':row[1]})
# Main loop
for feature in features:
gid = feature.attributes()[fieldIdx]
# if the feature has got an error
if gid in [x['GID'] for x in errors]:
error = (x['ERROR'] for x in errors if x['GID'] == gid).next()
geom = feature.geometry()
attrs = feature.attributes()
# write the feature to the output layer
outFeat.setGeometry(geom)
attrs.append(error)
outFeat.setAttributes(attrs)
writer.addFeature(outFeat)
current += 1
progress.setPercentage(int(current * total))
del writer
Register GeoAlgo
For QGis to display the GeoAlgo in the algorithms tree, you have to register the Python module first. Registering the module is done by editing the file named QGISAlgorithmProvider.py.
Just add the following lines:
...
# on the declaration of modules
from VerifyOracle import VerifyOracle
...
self.alglist = [SumLines(), PointsInPolygon(), # just add a line to this list
...
VerifyOracle(),
In action
Just open an Oracle layer in QGis, then launch the GeoAlgo. It is located in QGis GeoAlgorithms and in Vector analysis tools:
Then you'll find the following dialog box:
- In Vector layer choose the layer you want to verify.
- In UID Field for vector layer choose an attribute which is a primary key for the layer.
- In Result vector layer, you can add a path to a new shapefile or leave it blank to use a temporary file.
- Click on the Run button.
Once run and done, you can find a layer named "Result Vector Layer" that contains the non-valid geometries. Every attributes are taken from the original layer and one column has been added at the end, which is named Errorsand contains what errors have been found by Oracle. The ring and edges numbers are very useful to find where is the problem as QGis uses the same numbers in its nodes tool:
Conclusion
This code is provided as is ! There is plenty of room for improvement. You should consider it as something experimental... But it works !
Today I have been stuck by something quite strange under the Python environnement of QGis. I just wanted to open an Oracle vector layer and add it to the layer registry.
Here is what I tried (you can try it in the QGis Python Console from Extensions menu):
uri = QgsDataSourceURI()
uri.setConnection("","1521", "geobase", "user", "password")
uri.setDataSource("SCHEMA","TABLE",
v = QgsVectorLayer(uri.uri(), "MY_LAYER", "oracle")
QgsMapRegistry.instance().addMapLayer(v, True)
When launching the addMapLayer method, nothing moves on QGis interface and after a bit of digging in the log messages, I found nothing suspicious...
I tried to grab the state of the layer. You can do it with the isValid method which is inherited from the QgsMapLayer class. With all of the layers I tried to open, v.isValid() always returned False.
I tried to use the error() method which returned a QgsError object to try to find why my layer was invalid. But v.error().message() was always empty.
I remembered that I've got a special layer in my Oracle database. It is the only one which is well declared in Oracle metadata table (ALL_SDO_GEOM_METADATA view). When I tried the above lines code, I was able to open the layer in QGis map dialog. This special layer just have a well declared SRID, a well declared extent and stores only one type of geometry (Polygon).
I immediately struggle in QGis Python Console to find how to force a CRS for a QgsVectorLayer. The only way I found was to use the setCrs() method directly from the QgsVectorLayer class:
uri = QgsDataSourceURI()
uri.setConnection("","1521", "geobase", "user", "password")
uri.setDataSource("SCHEMA","TABLE",
v = QgsVectorLayer(uri.uri(), "MY_LAYER", "oracle")
v.setCrs(QgsCoordinateSystem(u"EPSG:27562")
But declaring the CRS is finally not required: if there is no SRID for the layer, QGis just use the CRS of the project (it depends on what you have configured in your QGis parameters). Furthermore, adding a CRS with the above method does not change the isValid() status...
Actually, the real problem was a wkbType one. All of my unopenable layers were declared with a wkbType to WkbUnknown. After re-reading QGis API documentation, the only way I've find so far to force the wkbType of a layer is to use uri.setWkbType() method. There is no other choice.
To open a layer with multiple geometry types, you have to add a little bit of code to the above snippet:
uri = QgsDataSourceURI()
uri.setConnection("","1521", "geobase", "user", "password")
uri.setDataSource("SCHEMA","TABLE",
uri.setWkbType(QGis.WkbPolygon)
v = QgsVectorLayer(uri.uri(), "MY_LAYER", "oracle")
QgsMapRegistry.instance().addMapLayer(v, True)
It took just 3 hours of work and search to find the solution. Why not make it public ?
Introduction
QGis can open a wide bunch of raster formats even proprietary ones. If you have ArcGis 9.3 (well, if you are forced to have it !), you have perhaps generated map services. Those services can be shown on cartographic web applications by using some kind of tile map service (TMS). This service from Esri ArcGIS doesn't respect the TMS standard. The TMS standard uses a special projection (EPSG:3857) and a well-known algorithm to determine which tile to load. Google Maps and OpenStreet map tiling services use the TMS standard and you are able to show them in QGis with just an XML description file (see this).
But it is not the case with ArcGIS 9.3 which uses a different algorithm to calculate which tile should be requested. This "no-standard" seems to have been discarded in the next version of ArcGIS (10.x). But after a bit of search, I finally manage to show an ArcGIS 9.3 map service on QGis 2.4 and that's what will be presented in this article.
Methodology
I've just studied what can open ArcGIS 9.3 cache maps and found that nearly nothing can do it, with the exception of ESRI ArcGIS of course ! Actually, in the opensource world, nothing else than OpenLayers 2 (starting to version 2.11) can open it ! I think that it is so exotic that it is not supported by OpenLayers 3. On the other side, ESRI seems to have changed its way of dealing with TMS in the version 10.x of ArcGIS. I haven't tested it myself but there are plenty of examples on the Internet that shows that QGis can open ArcGIS 10.x Map Services with only an XML description file (for those who want to go further, they can try the experimental "ArcGIS REST API connector" plugin).
Then I realised that there is a QGis plugin called "OpenLayers". It is used to show GoogleMaps/Bing/OSM tiles by using a workaround at the times where QGis (actually GDAL) was not able to handle those rasters. It is still installed in the version 2.4 of QGis even if the plugin has been refactored.
Here what I've found about "OpenLayers" QGis Plugin underwork:
- There is a core plugin (in Python) code.
- This code uses HTML description files which are a simple page which calls Javascript to print a Map. You can open those files in your Web browser and content will be shown.
- Core code executes Javascript and grabs the images.
- When you zoom or move onto the map, core code execute an Openlayers function and grabs the images.
Quite simple. It is true that this seems to be an ugly hack. Using a browser engine (thanks to Qt) to grab some raster files from the web is not what I call a web service. But it works !
In order to open your ArcGIS 9.3 Map Service, you have to write an HTML description file which uses OpenLayers 2 special functions for ArcGIS 9.3 and let do the magic of the plugin operates ! Actually it is a bit more complicated...
Add the good version of OpenLayers
First of all, you need the good version of OpenLayers 2. Grab the latest from the Internet and save it in the plugin repository: .qgis2/python/plugins/openlayers_plugin/weblayers/html. Name it OpenLayers2.js for disambiguation with other version of the library that are already installed in the directory (OpenLayers.js and OpenLayers-2.8.js). You need the full version of the Javascript code, don't use mobile or light version because they do not contain the ArcGIS code.
Which files should be modified ?
Before diving into the code, be sure to have a global vision of what to change or to add (the root of the plugin is .qgis2/python/plugins/openlayers_plugin):
- Modify 'weblayers/weblayer.py': Add a class to handle ArcGIS Cache services not projected in ESPG:3857.
- Create 'weblayers/html/layer_x.html': New HTML(+JS) file(s) to add access to one service from ArcGIS 9.3 Cache Map. Add a new file for every new service.
- Create 'weblayers/arcgis.py': New file to create a group of layers which will be named "ArcGIS Layers" in the plugin menu.
- Modify 'openlayers_plugin.py': Import the new group of layers in the plugin and add them to the menu.
Once everything has been made, your layers will be available under the menu Internet -> OpenLayers plugin.
Hack core code
Once you've got the good version of OpenLayers, you have to modify the core code of the plugin (but just a very little bit). The hack is just a workaround to deal with projection. If you inspect .qgis2/python/plugins/openlayers_plugin/weblayers/weblayer.py, you'll find a definition for the WebLayer class from which the code deals with the declaration of the layers. Particulary, there is code to handle the projection. Then you find a children class called WebLayer3857. This class is used in every python group declaration file (a file that just declares what are the name of your layers and which HTML file is linked). It just declares a layer with the projection EPSG:3857 which is the projection of all TMS services (read the introduction of this article for further details).
The problem with this static projection declaration is that the plugin cannot handle TMS that are not projected into this projection. Actually, the ArcGIS 9.3 services I used for this test QGis were using EPSG:27592 and even with the correct HTML definition map file, I was not able to show them in QGis if I use a WebLayer3857 declaration. The requested tiles were always too far away from what the service was able to provide.
As a workaround, we just have to declare a new class which is derived from WebLayer and which just says that the projection of the layers declared with this class is always 27592. Just add the following lines at the end of weblayer.py:
class ArcGISLayer27592(WebLayer):
epsgList = [27592]
The code of WebLayer class is just self-sufficient to handle the different projection. Then, you'll have to use the new ArcGISLayer27592 class to declare your layers in the group layer definition file.
Note to myself: the official plugin code for WebLayer could be updated to handle different kinds of projections...
Create the HTML definition map file
You just have to write down the definition of the map in an HTML file. Actually, it is much more Javascript than HTML. In the Javascript code, you'll have to use the OpenLayers functions that are able to open ArcGIS 9.3 Cache map. Be sure to load the good version of the OpenLayers2 library (OpenLayers2.js). Here is such a file:
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<title>ArcGIS Cache Map</title>
<link rel="stylesheet" href="qgis.css" type="text/css">
<script src="OpenLayers2.js"></script>
<script type="text/javascript">
var map;
var loadEnd;
var layerURL = "http://arcgis.example.com/arcgis/rest/services/folder_name/layer_name/MapServer";
var layerData = {
"serviceDescription": "",
"mapName": "Couches",
"description": "",
"copyrightText": "",
"layers": [
{
"id": 0,
"name": "The layer",
"parentLayerId": -1,
"defaultVisibility": true,
"subLayerIds": null
}
],
"spatialReference": {
"wkid": 27592
},
"singleFusedMapCache": true,
"tileInfo": {
"rows": 512,
"cols": 512,
"dpi": 96,
"format": "PNG32",
"compressionQuality": 0,
"origin": {
"x": 225000,
"y": 401000
},
"spatialReference": {
"wkid": 27592
},
"lods": [
{
"level": 0,
"resolution": 0.6614596562526459,
"scale": 2500
},
{
"level": 1,
"resolution": 0.26458386250105836,
"scale": 1000
},
{
"level": 2,
"resolution": 0.13229193125052918,
"scale": 500
},
{
"level": 3,
"resolution": 0.06614596562526459,
"scale": 250
}
]
},
"initialExtent": {
"xmin": 301177.13180509704,
"ymin": 258249.61360172715,
"xmax": 301317.7581280163,
"ymax": 258352.66901617133,
"spatialReference": {
"wkid": 27592
}
},
"fullExtent": {
"xmin": 270184,
"ymin": 235952,
"xmax": 332136,
"ymax": 281008,
"spatialReference": {
"wkid": 27592
}
},
"units": "esriMeters",
"supportedImageFormatTypes": "PNG32,PNG24,PNG,JPG,DIB,TIFF,EMF,PS,PDF,GIF,SVG,SVGZ",
"documentInfo": {
"Title": "Layer",
"Author": "author",
"Comments": "",
"Subject": "",
"Category": "",
"Keywords": "",
"AntialiasingMode": "Normal",
"TextAntialiasingMode": "Force"
}
};
function init() {
loadEnd = false;
function layerLoadStart(event)
{
loadEnd = false;
}
function layerLoadEnd(event)
{
loadEnd = true;
}
var baseLayer = new OpenLayers.Layer.ArcGISCache(
"AGSCache",
layerURL,
{
layerInfo: layerData,
eventListeners: {
"loadstart": layerLoadStart,
"loadend": layerLoadEnd
}
}
);
map = new OpenLayers.Map('map', {
maxExtent: baseLayer.maxExtent,
units: baseLayer.units,
resolutions: baseLayer.resolutions,
numZoomLevels: baseLayer.numZoomLevels,
tileSize: baseLayer.tileSize,
displayProjection: baseLayer.displayProjection
});
map.addLayer(baseLayer);
}
</script>
</head>
<body onload="init()">
<div id="map"></div>
</body>
</html>
You can see that there is a big variable called layerData. It describes the whole content definition of the ArcGIS 9.3 service. Without this variable, you are not able to open the layer in an OpenLayers map... you can grab the variable value by using this example from OpenLayers2 which uses the autoconfig mechanism of ArcGISCache. Just open the Firefox Developper Tools and get the variable value (or use console.info()).
Howto grab layer definition easily ?
Creating HTML definition map file is something that can be very fastidious. Actually, you have to grab the JSON value for the layerData variable which describes the content of the service. You can do it by using Firefox Developper Tools but this is not an effective method. To ease the pain, just load this HTML (+JS) file in your web browser and put the URL of an ArcGIS 9.3 Map service (for example: http://arcgis.example.com/arcgis/rest/services/name_of_folder/name_of_layer/MapServer) before clicking on the load button. You'll have to save the resulting HTML file in the plugin directory: .qgis2/python/plugins/openlayers_plugin/weblayers/html.
You need to have Internet access to use this code because I just grab OpenLayers2 from cloudflare CDN. But you can replace it with any other javascript file that contains OpenLayers.
<!DOCTYPE html>
<html>
<head>
<title>QGis OpenLayers ArcGIS 9.3 Cache config helper</title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<style>
.smallmap {
width: 712px;
height: 556px;
border: 1px solid #ccc;
}
</style>
<script src="http://cdnjs.cloudflare.com/ajax/libs/openlayers/2.11/OpenLayers.js"></script>
<script type="text/javascript">
var map;
var layerURL;
function init() {
var url = document.getElementById("url").value;
if (url) {
layerURL = url;
var jsonp = new OpenLayers.Protocol.Script();
jsonp.createRequest(layerURL, {
f: 'json',
pretty: 'true'
}, initMap);
}
}
function initMap(layerInfo){
var layerData = JSON.stringify(layerInfo, null,'\t');
var html_doc = "<html xmlns=\"http://www.w3.org/1999/xhtml\">\n\
<head>\n\
<title>ArcGIS Cache Map<\/title>\n\
<link rel=\"stylesheet\" href=\"qgis.css\" type=\"text\/css\">\n\
<script src=\"OpenLayers2.js\"><\/script>\n\
<script type=\"text/javascript\">\n\
var map;\n\
var loadEnd;\n\
var layerURL = \"" + layerURL + "\";\n\
var layerData = " + layerData +";";
html_doc = html_doc + "\n\
function init() {\n\n\
loadEnd = false;\n\
function layerLoadStart(event)\n\
{\n\
loadEnd = false;\n\
}\n\
function layerLoadEnd(event)\n\
{\n\
loadEnd = true;\n\
}\n\
var baseLayer = new OpenLayers.Layer.ArcGISCache(\n\
\"AGSCache\",\n\
layerURL,\n\
{\n\
layerInfo: layerData,\n\
eventListeners: {\n\
\"loadstart\": layerLoadStart,\n\
\"loadend\": layerLoadEnd\n\
}\n\
}\n\
);\n\n\
map = new OpenLayers.Map('map', {\n\
maxExtent: baseLayer.maxExtent,\n\
units: baseLayer.units,\n\
resolutions: baseLayer.resolutions,\n\
numZoomLevels: baseLayer.numZoomLevels,\n\
tileSize: baseLayer.tileSize,\n\
displayProjection: baseLayer.displayProjection\n\
});\n\
map.addLayer(baseLayer);\n\
}\n\
<\/script>\n\
<\/head>\n\
<body onload=\"init()\">\n\
<div id=\"map\"><\/div>\n\
<\/body>\n\
<\/html>";
var html_textarea = document.getElementById("qgis_html");
html_textarea.value = html_doc;
// Send the HTML file to save in the local filesystem:
var a = document.createElement('a');
a.href = 'data:attachment/html,' + encodeURIComponent(html_doc);
a.target = '_blank';
a.download = 'arcgis_layer.html';
document.body.appendChild(a);
// send it to the browser
a.click();
var baseLayer = new OpenLayers.Layer.ArcGISCache("AGSCache", layerURL, {
layerInfo: layerInfo
});
/*
* Make sure our baselayer and our map are synced up
*/
map = new OpenLayers.Map('map', {
maxExtent: baseLayer.maxExtent,
units: baseLayer.units,
resolutions: baseLayer.resolutions,
numZoomLevels: baseLayer.numZoomLevels,
tileSize: baseLayer.tileSize,
displayProjection: baseLayer.displayProjection
});
map.addLayer(baseLayer);
map.addControl(new OpenLayers.Control.LayerSwitcher());
map.addControl(new OpenLayers.Control.MousePosition() );
var bbox = layerInfo.initialExtent;
map.zoomToExtent(new OpenLayers.Bounds(bbox.xmin, bbox.ymin, bbox.xmax, bbox.ymax));
}
</script>
</head>
<body>
<h1 id="title">OpenLayers ArcGIS 9.3 Cache config helper</h1>
<p id="shortdesc">
<form id="ServiceURL">
<p>
<label for="url">Service URL:</label>
<input id="url" type="url" style="width: 30%;" placeholder="http://arcgis.example.com/arcgis/rest/services/folder_name/layer_name/MapServer"/>
<input type="button" value="Load" onclick="init();"/>
</p>
</form>
</p>
<form id="serviceFolders">
</form>
<div id="map" class="smallmap"></div>
<div id="docs">
<p>HTML to save in a file</p>
<textarea id="qgis_html" style="width: 40%; height: 300px;"></textarea>
</div>
</body>
</html>
Make a group definition file
You have to build a group definition file. This is a purely declarative python file which just describes what are your layers. It is used to build the menu shown by the plugin and to make the link with the HTML files. The first class is the Group one. Then all of the layers that are described in HTML file have to be children class of the group one. Be sure to use strings (and not unicode things) in the name of the layers.
# -*- coding: utf-8 -*-
from weblayer import ArcGISLayer27592
class OlArcgisLayer(ArcGISLayer27592):
emitsLoadEnd = True
def __init__(self, name, html):
ArcGISLayer27592.__init__(self, groupName="ArcGIS 9.3", groupIcon="arcgis_icon.png",
name=name, html=html)
class OlArcgisFondDePlanLayer(OlArcgisLayer):
def __init__(self):
OlArcgisLayer.__init__(self, name="Fond de plan", html="arcgis_fond_de_plan.html")
class OlArcgisEquipementPublicLayer(OlArcgisLayer):
def __init__(self):
OlArcgisLayer.__init__(self, name='Equipement public', html='arcgis_equipement_public.html')
class OlArcgisAdresseLayer(OlArcgisLayer):
def __init__(self):
OlArcgisLayer.__init__(self, name='Adresse', html='arcgis_adresse.html')
class OlArcgisArreteCirculationLayer(OlArcgisLayer):
def __init__(self):
OlArcgisLayer.__init__(self, name='Arrete de circulation', html='arcgis_arrete_circulation.html')
class OlArcgisAssainissementLayer(OlArcgisLayer):
def __init__(self):
OlArcgisLayer.__init__(self, name='Assainissement', html='arcgis_assainissement.html')
class OlArcgisCadastreNapoleonLayer(OlArcgisLayer):
def __init__(self):
OlArcgisLayer.__init__(self, name='Cadastre Napoleonien', html='arcgis_cadastre_napoleon.html')
class OlArcgisParcelleLayer(OlArcgisLayer):
def __init__(self):
OlArcgisLayer.__init__(self, name='Parcelle', html='arcgis_parcelle.html')
Register group definition file in the plugin
You have to modify openlayers_plugin.py to add the entry menus from your group definition file. The entries are built upon the python class, so you have to import them with a simple line like:
...
from weblayers.arcgis import OlArcgisFondDePlan, OlArcgisEquipementPublic, ...
...
Once the modules (classes) are declared, you need to manually add them into the function initGui:
...
self._olLayerTypeRegistry.register(OlArcgisFondDePlanLayer())
self._olLayerTypeRegistry.register(OlArcgisEquipementPublicLayer())
...
Once this is done, you should be able to find your menu entries in the menu with a group of your own layers. Be careful about the names which can be a little type error prone.
User experiment
The OpenLayers QGis plugin is far from being perfect ! Actually it seems to be an ugly hack that works with very few efforts, so on... I've experienced some troubles by using this modified plugin with ArcGIS 9.3:
- You should open another layer which is in the extent of your raster layer before opening the raster layer. Otherwise, the HTML script may crash because you are asking for a too much wide extent. It can block QGis, so be careful !
- Sometimes when you zoom, not of all the tiles are shown on the map dialog. You have to unzoom and rezoom to have the whole map updated and shown.
- Same problem with panning or moving. There you have to unzoom and rezoom to grab all of the tiles.
- It seems to have a sort of cache: once you grabbed all the tiles at different zoom levels, map loading is more fast.
- You can't have more than one OpenLayers layer on the map. Only the latest opened layer is shown, whatever the layer order. It is a limitation of the whole plugin, not restricted to ArcGIS services.
- You can't save an OpenLayer layer on your project ! You have to reopen the layer each time you open your project. It is also a limitation of the plugin, not restricted to ArcGIS services.
I have taken time to test print composer with those services and it seems to be unaffected by the bugs I've seen in the map window. You will be able to generate PDF or print maps with your ArcGIS Map service layers. So far, so good...
Conclusion
QGis is good for you. It can even be interfaced with its greatest opponent in GIS Desktop: ESRI ArcGIS. If you haven't migrated to ArcGIS 10.x, you are still able to use your old ArcGIS 9.3 Cache map services for a while under QGis. There are great chances that no other GIS Desktop software is able to deal with these services nowadays. But QGis is one of them... To go further, I can only advice you to migrate to a true free software tiling solution like MapNik and carto-css. By using the TMS standard, you are 100% sure that QGis will be able to do it. Please read my wiki article to build your opinion on the ease to generate tiles with those tools.