Représenter un résultat de requêtes spatiales sous PostgreSQL PostGIS depuis un client web

Dans l’objectif de développer une application cartographique, nous avons au niveau des données récupéré des jeux de données spatiales, installé la base de données et son extension spatiale et finalement chargé les données en base. Nous avons également mis en place un serveur cartographique qui fournit les données au travers du réseau, et construit les bases de l’application dans un client web pour interroger le serveur cartographique et afficher les premiers éléments de carte.

L’application ne se contente pas d’afficher des données statiques, mais elle doit permettre de représenter l’information dynamiquement en effectuant des requêtes en fonction du contexte d’utilisation. C’est la partie que nous allons développer dans ce document.

Préparation de l’interface cartographique

La première étape est de préparer l’interface graphique et recadrer le contexte.

Nous avons vu comment appeler les services WMS et WFS que nous avons mis en place pour générer une carte simple où la proximité des établissements scolaires aux parcelles cultivées sur la Gironde était affichée de façon statique sur une carte rendue disponible dans un navigateur web par les librairies Leaflet et OpenLayers 3.

Dans la suite de ce document j’utiliserai la librairie Leaflet mais le principe est identique pour OpenLayers 3 et le code peut être adapté.

Question de volumétrie nous n’avons pas sur ce prototype intégré les données parcellaires sur toute les régions. Voyons en premier lieu comment adapter le code pour afficher les parcelles cultivées sur la France métropolitaine, puis explorons les fonctionnalités à paramétrer sur la carte.

Représentations des couches cartographiques depuis le serveur cartographique

Premièrement le programme représente le fond de carte par l’appel aux serveurs de tuiles d’OpenStreetMap, puis il fait appel au serveur QGIS Server pour récupérer les données géométriques des départements français dont il effectue le rendu à la volée. Ces fonctionnalités sont identiques à celles du premier prototype.

Ensuite le programme effectue le rendu des parcelles pour la France métropolitaine. Il interroge le serveur WMS qui génère une image pour la vue en cours.

var wmsLayerParcelles = L.tileLayer.wms(project, {
    layers: 'Parcelles',
    format: 'image/png',
    transparent: true,
    minZoom: 10,
    maxZoom: 18,
    opacity: 0.55,
});
wmsLayerParcelles.addTo(map);

Cette fonctionnalité a été expliquée, la nouveauté est l’apparition des paramètres minZoom et maxZoom qui vont permettre de conditionner le rendu de la couche au niveau d’agrandissement en cours sur la carte. En d’autres termes la librairie ne fera pas appel au serveur pour représenter les parcelles en dessous d’un niveau d’agrandissement fixé à 10. Cela permet de limiter la quantité d’informations de parcelles sur une image, et donc à la fois le travail de représentation des données du serveur et le trafic sur le réseau où seraient transférées des images trop volumineuses.

Si le rendu est effectué côté serveur, on souhaite également accéder aux caractéristiques des parcelles pour indiquer le numéro d’îlot et le type de culture. Cela peut être fait à un niveau d’agrandissement encore supérieur afin de limiter la taille des fichiers échangés, en faisant appel au service WFS.

map.on('zoomend moveend', function() {
...
    if (map.getZoom() > 15)
        getJsonLayerParcelles(map.getBounds());
});

On va préciser au gestionnaire d’événements de mettre en écoute les opérations d’agrandissement et de glisser-déposer.

Lorsque ces événements arrivent, on vérifie le niveau d’agrandissement de la carte et s’il est supérieur à la valeur seuil que l’on a choisi pour faire figurer les données, on charge la couche vectorielle de façon asynchrone pour l’étendue géographique en cours sur la carte.

var jsonLayerParcelles = null;

function getJsonLayerParcelles(bounds) {

    var getFeatureParcelles = project + '?SERVICE=WFS'
        + '&VERSION=1.1.0'
        + '&REQUEST=GetFeature'
        + '&TYPENAME=Parcelles'
        + '&SRSNAME=EPSG:4326'
        + '&outputFormat=GeoJSON'
        + '&BBOX=' + bounds.toBBoxString();

    var proxy = host + '/proxy.php?callback=getGeoJson&url=';

    var proxyURL = proxy + encodeURIComponent(getFeatureParcelles);

    var jsonLayerParcelles = L.geoJson(null, {
        style: function (feature) {
            return {
            fillColor: 'rgb(255,255,255)',
            color: 'rgb(175,179,138)',
            weight: 1,
            opacity: 1,
            fillOpacity: 0,
            dashArray: '1,5',
            };
        },
        onEachFeature: function (feature, layer) {
            var content = 'Parcelle N° ' + feature.properties['num_ilot'] + '<br />';
            content += getCulture(feature.properties['cult_maj']);
            layer.bindPopup(content);
        }
    });

    $.ajax({
        url: proxyURL,
        dataType: 'jsonp',
        jsonpCallback: 'getGeoJson',
        success: function (response) {
            L.geoJson(response, {
                onEachFeature: function (feature, layer) {
                    jsonLayerParcelles.addData(feature)
                }
            });
        }
    });
    jsonLayerParcelles.addTo(map);
}

Le principe de chargement de la couche vectorielle est le même que pour les départements sauf que cette fois-ci on passe l’étendue géographique dans l’url du service. Cela se fait via deux opérations, map.getBounds() et bounds.toBBoxString().

Remarquez que pour chaque objet récupéré on attache un contenu qui s’affichera dans un popup au clic sur la couche. Le contenu est mis en forme par rapports aux propriétés récupérées.

Comme le rendu est déjà effectué sous forme d’image côté serveur, l’opacité est fixée pour rendre transparent le remplissage des polygones, le contour est quand à lui représenté en pointillés grâce à l’option de style dashArray. Cela permet de faire figurer les parcelles qui ont été représentées en blanc opaque ou transparence car ce ne sont pas des cultures mais des prairies ou des gels de cultures.

Finalement par rapport au premier prototype, un contrôle a été ajouté pour utiliser le géocodage d’adresses afin de localiser une destination sur la carte.

        <link rel="stylesheet" href="css/Control.OSMGeocoder.css" />
        <script src="js/Control.OSMGeocoder.js"></script>
var osmGeocoder = new L.Control.OSMGeocoder({
    collapsed: false,
    position: 'topright',
    text: 'Allez!'
});
osmGeocoder.addTo(map);

Paramétrage de fonctionnalités sur la carte

Une extension à la librairie Leaflet permet de faire figurer sur la carte une barre latérale qui a la capacité de se déplier et replier au clic sur les icônes de la barre.

J’ai adapté les styles css pour effectuer un rendu qui laisse entrevoir le fond de carte par transparence.

L’extension est disponible sous Leaflet comme sous Openlayers, elle me semble un bon choix relativement standardisé pour apporter des fonctionnalités à la carte sans développer une multitude de contrôles.

L’extension requiert le chargement de ressources css et javascript, ainsi que du jeux d’icônes vectorielles sous police de caractères font-awesome. Je fais également appel aux polices de caractères Google pour le rendu du contenu de la barre latérale.

        <link href='http://fonts.googleapis.com/css?family=Lato:400,700|PT+Sans:400,700|Roboto:400,500,700' rel='stylesheet' type='text/css'>
        <link href="http://maxcdn.bootstrapcdn.com/font-awesome/4.1.0/css/font-awesome.min.css" rel="stylesheet">
...
        <script src="js/leaflet-sidebar.min.js"></script>

L’extension requiert de modifier la classe de style appliquée au conteneur de carte.

<div id="map" class="sidebar-map"></div>

Le contenu de la barre est simplement mis en place par une convention de balises html que je ne reprendrai pas ici, l’extension étant bien documentée.

Un onglet de paramètres permettra à l’utilisateur de configurer la carte. En termes de fonctionnalités on souhaite à minima paramétrer la distance qui caractérise la proximité des établissements aux parcelles, d’autres fonctionnalités peuvent permettre de faire une sélection des établissements à représenter en fonction de la surface des parcelles en contact, des types de cultures, etc…

Ces fonctionnalités requierent d’interroger la base de données après soumission d’un formulaire, l’interrogation se fait au moyen de requêtes spatiales sur Postgis, sachant que l’appel s’effectue en Ajax et que le serveur retournera des informations en Json.

Requêtes spatiales

Avant de passer au mécanisme de communication entre la carte et la base de données, nous allons procéder à l’étude des requêtes spatiales qui feront toute l’importance des données représentées.

La détermination de risques est une problématique complexe, ce que l’application propose n’est pas d’évaluer un risque mais d’exposer des facteurs de risques potentiels qui découlent simplement du bon sens. La démarche est de déterminer des populations d’élèves qui pourraient être impactées par ces critères afin de construire des échantillons de population sur lesquels effectuer des prélèvements pour évaluer l’exposition effective aux pesticides.

Intuitivement le risque d’exposition aux pesticides augmente quand :

  • L’aire cultivée augmente à proximité de la population. Une requête spatiale doit opérer un calcul d’aire. L’aire peut être représentée par un symbole proportionnel.
  • Il existe plusieurs types de cultures différentes à proximité d’une population. (multiplication de la nature des pesticides, et des horaires d’épandage). Une requête doit permettre de comptabiliser les cultures à proximité d’une population
  • La catégorie de culture est plus ou moins consommatrice de pesticides. Un indicateur existe, l’IFT, qui permet de mesurer la pression phytosanitaire (quantités et fréquences d’épandage) d’un sol mais il ne prend pas en compte la toxicité des pesticides, le type de culture sera par conséquent affiché seulement à titre informatif, il donne quand même une information exploitable car on a connaissance par exemple de la toxicité des traitements en viticulture. Les requêtes doivent éliminer les pâturages et gels de cultures

La carte sera imparfaite. Un certain nombre de paramètres influeront les résultats, par exemple :

  • Il n’est pas possible de prendre en compte la rotation des cultures. Les données de parcelles sont des déclarations relevées en 2010, ces données sont anciennes, et comportent nécessairement des erreurs. Un parallèle pourra être mis en place avec une carte aérienne.
  • Les cultures bio sont négligées faute de données. (De l’ordre de 3% des cultures)
  • Le centroide des établissements n’est pas l’aire des bâtiments, si l’on considère une proximité de 500 mètres aux coordonnées de latitudes et longitude de l’établissement, la proximité réelle aux salles de classe ou à la cours de récréation peut être moindre
  • Les vents dominants et l’age du proviseur ne sont pas pris en compte !

Fonctions spatiales PostGIS, géométrie ou géographie ?

Consultez la référence PostGIS !

PostGIS dispose d’un jeu très complet de fonctions spatiales, et l’on abordera que très peu d’entre elles. Sachez toutefois que PostGIS change de convention de nommage pour standardiser le nom des fonctions par l’ajout du préfixe spatial ST_ à chaque fonction. Cela signifie que si vous recherchez des exemples d’utilisation sur Internet il vous faudra probablement retirer le préfixe…

Les fonctions spatiales peuvent pour certaines fonctionner avec plusieurs types de paramètres. Des conversions (implicites ou explicites) peuvent être effectuées entre certains types. Les types de paramètres sont :

  • box2d : est une boîte composée de coordonnées xmin, ymin, xmax, ymax. Souvent utilisé pour retourner la boîte englobante 2D d’une géométrie.
  • box3d : est une boîte composée de coordonnées xmin, ymin, zmin, xmax, ymax, zmax. Souvent utilisé pour retourner la mesure 3D d’une géométrie ou le recouvrement des géométries.
  • geometry : est le type de données spatiales utilisé par PostGIS pour représenter une fonction dans le système de coordonnées Euclidien. L’unité est celle du SRID.
  • geometry_dump : est un type de données spatiales avec deux champs, geom et path[] qui sont respectivement un objet de géométrie et un tableau qui contient la position de la géométrie dans l’objet de dump
  • geography : est un type de données spatiales utilisée pour représenter une entité dans le système de coordonnées sphériques terrestres. L’unité est le mètre carré.

Nous utiliserons les types geometry et geography. Nos données sont enregistrées dans le type geometry et pour nos besoins nous voudrons calculer des distances en mètres ou des superficies en mètres carrés, cela implique de faire des conversions de type geometry vers type geography et vice versa.

geom::geography
geom::geometry

Sont des exemples de conversions explicites. (Vous pouvez vous référer à la documentation pour savoir quels types peuvent être convertis en autres types)

Dernière remarque concernant les types, l’unité géographique d’une superficie est le mètre carré. Pour des raisons évidentes nous emploieront plus volontiers le kilomètre carré ou l’hectare pour donner une information.

Les conversions d’unités sont les suivantes :

1 mètre carré = 1 * 10^-6 kilomètre carré   (ou encore 1 / 1000000)
1 mètre carré = 1 * 10^-4 hectare   (ou encore 1 / 10000)

Préparation des requêtes

Calcul de l’aire d’un polygone

ST_Area – Retourne l’aire de la surface d’un polygone ou multi-polygone.

Pour l’exemple je vais calculer l’aire de chaque département. La colonne est de type géométrique, il suffit de faire une conversion de type pour retourner la surface en mètres carrés.

SELECT nom_dept, ST_Area(CAST(geom As geography)) * POWER(10, -6) as geom FROM "Departement";

Distance entre géométries

La recherche d’établissements à proximité de zones cultivées peut être réalisée au travers de plusieurs fonctions, le problème qui va se poser est celui des performances.

ST_DWithin – Renvoie la valeur vraie si les géométries sont dans la distance spécifiée l’une à l’autre. Pour les unités géographiques la distance est exprimée en mètres, la mesure s’effectue par défaut autour d’un sphéroïde.

ST_Distance – Pour des champs de type géographie, retourne la distance minimale entre deux géographies sphéroïdale en mètres.

ST_Buffer – Pour un champ de type géométrie, retourne une géométrie qui représente tous les points dont la distance à l’objet géométrique est inférieure ou égale à la distance. Les calculs sont dans le système de référence de l’objet. Pour une géométrie de type point le résultat sera approximativement un cercle ayant pour centre le point et pour rayon la distance.

ST_Intersects – Retourne la valeur vraie si les géométries / géographie « se croisent spatialement en 2D » – (possèdent une partie de l’espace en commun) et la valeur faux sinon (les objets sont disjoints). Pour la géographie – la tolérance est 0,00001 mètres (donc tous les points qui sont considérés fermés retournent la valeur vraie)

Si l’on s’appuie sur les deux premières fonctions, il suffit en théorie de faire une jointure de la table d’établissement avec la table (ou la vue matérialisée) des cultures sur la géométrie des objets trans-typée en géographie sur la distance souhaitée. L’approche est intuitive mais n’exécutez pas la requête qui suit…, pour une proximité de 50 mètres on peut tenter d’évaluer le temps d’exécution :

EXPLAIN SELECT 
e.numero_uai, e.denominati, e.geom, c.num_ilot, c.cult_maj
FROM "Etablissement" e, "master_cultures" c 
WHERE ST_DWithin(e.geom::geography, c.geom::geography, 50);

Le coût d’exécution est estimé à :

"Nested Loop  (cost=0.00..119895256281.96 rows=132693 width=104)"

La requête prendrait un temps indécent à s’exécuter !

Pour une utilisation en ligne il est impératif de compiler au préalable les données. Ma première idée est d’établir un distancier entre le centroide des établissements et les parcelles de cultures à moins de 1500 mètres, le seuil de distance maximale qui sera proposé à l’utilisateur.

En volumétrie, nous avons 64.901 établissements et 6.132.686 parcelles, soit sans seuil potentiellement 398.017.454.086 lignes !

EXPLAIN 
SELECT count(e.numero_uai)
FROM "Etablissement" e, "master_cultures" c
WHERE ST_DWithin(e.geom::geography, c.geom::geography, 1500);
"Aggregate  (cost=119116282174.69..119116282174.70 rows=1 width=9)"

Cette approche permettrait d’enregistrer les distances dans une table indexée mais cette opération est toujours indécente si l’on considère le temps d’exécution initial.

Une telle quantité de parcelles cultivées est très pénalisante, peut-être peut-on éliminer certaines parcelles de l’équation ? Le type de culture est donné par un numéro. On peut éventuellement éliminer certains types, comme les prairies, qui ne sont pas à priori à prendre en considération.

select count(num_ilot) from master_cultures
where cult_maj not in (18, 19);

Retirer les prairies et prairies temporaires permet de diviser par deux le nombre de lignes de la table de cultures, cela reste toutefois très insuffisant.

Il faut envisager une autre approche. Si le coût d’une opération ST_DWithin() est trop important peut être que l’opération ST_Intersects() se révélera plus rapide.

Si l’on se base sur le centroide des établissements pour créer une zone tampon de 1500 mètres autour de l’établissement on peut chercher l’intersection entre la zone tampon et les parcelles, puis on pourra rechercher la distance sur les cas limités où la proximité est avérée.

create table tampon_etas as
select e.numero_uai, e.geom as centre_geom, ST_Buffer(e.geom::geography, 1500)::geometry as geom 
FROM "Etablissement" e;

La première opération est triviale, il suffit d’employer la fonction ST_Buffer() autour de la géométrie de l’établissement, de type POINT, sur un rayon de 1500 mètres pour créer une nouvelle géométrie. Cette opération prend environ 13 secondes, suite à quoi on peut ajouter des indexes spatiaux. La géométrie du point est conservée, elle matérialise le centre de la zone tampon, à partir duquel sera calculée la distance aux cultures.

CREATE INDEX tampon_cgeom_gist ON tampon_etas USING gist (centre_geom) TABLESPACE pg_default;
CREATE INDEX tampon_geom_gist ON tampon_etas USING gist (geom) TABLESPACE pg_default;

Tentons d’évaluer l’opération d’intersection…

EXPLAIN SELECT t.numero_uai
FROM "tampon_etas" t, "master_cultures" c 
WHERE ST_Intersects(t.geom, c.geom);
"Nested Loop  (cost=0.28..176252434.15 rows=132693478 width=9)"

Le coût est très élevé mais déjà divisé par 1000 par rapport aux premières approches. Nous pouvons tenter d’exécuter une requête sur les cultures d’un seul département, et multiplier le temps d’exécution par le nombre de départements pour avoir une estimation du temps d’exécution.

SELECT t.numero_uai
FROM "tampon_etas" t, "cultures_01" c 
WHERE ST_Intersects(t.geom, c.geom)
AND cult_maj not in (18, 19);

Malheureusement l’opération retourne une erreur !

ERREUR:  GEOSIntersects: TopologyException: side location conflict (...)

L’opération n’abouti pas car certaines données de la table de cultures sont invalides, PostGIS possède des fonctions qui permettent de vérifier les données, ou encore de réparer des données.
Si le taux d’erreur n’est pas trop élevé on pourra ignorer les parcelles dont les données sont en erreur.

ST_IsValid – Teste une valeur géométrique et retourne la valeur vrai si la géométrie est bien formée.

ST_IsSimple – Teste une valeur géométrique et retourne la valeur vrai si la géométrie n’a pas de point géométrique invalide comme une auto-intersection ou auto-tangence.

Notez que la géométrie de la table des cultures est de type MULTIPOLYGON, mais que certaines de ses entrées sont de type GeometryCollection, il y a une incohérence de type pour certaines entrées.

SELECT count(*)
FROM "cultures_01" c 
WHERE cult_maj not in (18, 19)
AND NOT ST_IsValid(geom);

142 géométries se révèlent invalides sur 63962 géométries testées, soit environ 2 pour mille. Par rapport à notre besoin on peut se contenter d’ignorer les parcelles dont la géométrie est invalide, le résultat sera négligeable.

SELECT t.numero_uai, c.num_ilot,
t.geom as eta_geom, c.geom as cul_geom,
ST_Distance(t.centre_geom::geography, c.geom::geography) as distance
FROM "tampon_etas" t, "master_cultures" c
WHERE
ST_IsValid(c.geom)
AND ST_Intersects(t.geom, c.geom)
LIMIT 100;

Lorsque l’intersection entre la zone tampon et les cultures retourne la valeur vrai on peut récupérer la distance du centre du tampon à la géométrie de la culture. Cela souffre un peu d’imprécision, mais la fonction ST_DWithin() est contre performante dans ce cas. Il ne reste plus qu’à créer la table qui conservera les distances, c’est cette table que nous interrogerons.

explain CREATE TABLE intersect_etabs AS 
SELECT t.numero_uai, e.code_dept, c.num_ilot,
t.centre_geom as eta_tampon_geom, c.geom as cul_geom,
CAST(floor(ST_Distance(t.centre_geom::geography, c.geom::geography)) as smallint) as distance
FROM 
	"tampon_etas" t INNER JOIN "Etablissement" e on t.numero_uai = e.numero_uai, 
	"master_cultures" c
WHERE
c.cult_maj not in (18, 19) 
AND ST_IsValid(c.geom)
AND ST_Intersects(t.geom, c.geom);

La création prend 25 minutes… J’ai transformé la distance calculée afin d’utiliser un entier de petite dimension qui sera plus rapide à interroger qu’un nombre avec décimales.

Une jointure sur la table d’établissements permet de récupérer au passage le code du département.

select * from intersect_etabs where distance < 300;

La requête prend environ 33 secondes sur l’ensemble de la table sans utiliser d’étendue. Ce résultat n’est pas exploitable sur l’étendue de la France métropolitaine, d’une part parce que le résultat est trop long à retourner, d’autre part parce que cela ferait représenter trop d’établissements sur la carte.

La solution est d’exécuter la requête uniquement lorsque le niveau d’agrandissement est adéquat, on profite alors des capacités de l’indexe spatial. Au niveau de la vue sur la France entière, on pourra faire afficher au niveau de chaque département le nombre d’établissements à proximité de parcelles cultivées dans un symbole de cercle proportionnel. Nous devons au préalable compiler les informations pour les départements pour chaque option de distance proposée à l’utilisateur dans l’interface.

Vérifions en premier lieu la réactivité de la requête spatiale sur une étendue géographique limitée.

Limiter la requête spatiale PostGIS à une étendue géographique

ST_MakeEnvelope – Crée un polygone rectangulaire formé à partir des minimums et maximums donnés. Les valeurs en entrée doivent être spécifiées dans le système de coordonnées géographiques spécifié par le SRID.

Sous PostGIS on peut utiliser la fonction ST_MakeEnvelope(left, bottom, right, top, srid) pour construire l’enveloppe de sélection de l’étendue géographique, qui associée à l’opérateur && sur la géométrie de l’objet permet de trouver l’intersection entre l’enveloppe et la géométrie.

CREATE INDEX intersect_etabs_cul_geom_gist
  ON intersect_etabs
  USING gist
  (cul_geom);
SELECT * FROM intersect_etabs 
WHERE cul_geom && ST_MakeEnvelope(0.23036956787109372,44.9590891448628,0.3277873992919922,45.00419734261587, 4326)
AND distance < 1000;

La requête s’exécute en 11 ms, ce qui est compatible avec une utilisation en temps réel sur la carte.

Compilation des données par département

L’objectif est de créer une table, relativement longue à générer mais très rapide à interroger.

WITH   proximites(a) AS ( VALUES ('{100, 250, 500, 750, 1000, 1250, 1500}'::int[]) )
SELECT generate_subscripts(a, 1) AS idx, unnest(a) AS proximite
FROM   proximites;

Cette première requête permet de générer l’ensemble des proximités aux établissements proposées à l’utilisateur, elle peut être utilisée dans une requête plus complexe de création de table.

CREATE TABLE proximites AS (
WITH   proximites AS (
WITH generateur(a) AS ( VALUES ('{100, 250, 500, 750, 1000, 1250, 1500}'::int[]) )
SELECT generate_subscripts(a, 1) AS idx, unnest(a) AS valeur
FROM   generateur
)
SELECT p.valeur as proximite, count(distinct i.numero_uai) as intersections, i.code_dept 
FROM proximites p
LEFT JOIN intersect_etabs i
ON i.distance < p.valeur
GROUP BY p.valeur, i.code_dept
ORDER BY p.valeur, i.code_dept
)

L’évolution de la table pour d’autres valeurs est simple à réaliser…

INSERT INTO proximites (
SELECT '250'::text as proximite, count(distinct i.numero_uai) as intersections, i.code_dept 
FROM intersect_etabs i 
WHERE i.distance < 250 
GROUP BY i.code_dept
)

Vérifions que les résultats sont conformes aux attentes de réactivité.

select intersections, code_dept from proximites where proximite = '250';

Nous sommes bien sous les 10 ms.

Nous avons compilé les données et examiné les requêtes spatiales qui vont permettre l’interactivité avec l’utilisateur depuis la carte. Nous pouvons désormais utiliser la table des intersections entre établissements et cultures sur une étendue géographique donnée, et à niveau de détail moins élevé la table qui recense le nombre d’établissements en fonction de la proximité aux cultures pour chaque départements.

Voyons comment représenter dynamiquement ces informations sur la carte.

Représenter les résultats JSON de requêtes spatiales

Le principe est simple. Une carte est crée qui constitue la vue par défaut. Une interaction utilisateur va déclencher le changement de la vue pour s’adapter au contexte. Techniquement un écouteur d’événements est mis en place dans le code Javascript qui s’exécute côté client sur le navigateur Internet.

Au déclenchement d’un événement, un appel Javascript asynchrone interroge un script distant qui retourne l’information au format JSON. Le script distant a pour responsabilité de faire des requêtes spatiales sur la base de données et mettre en forme le contenu au format d’échange.

Le contenu retourné peut comporter des informations spatiales à représenter sous la librairie Javascript, dans ce cas le contenu d’échange sera du JSON, qui pourra inclure des objets GeoJSON. Nous utiliserons ce principe pour récupérer les points qui permettrons de localiser et représenter les établissements dans la vue détaillée, toutefois ce n’est pas obligatoire : en effet dans notre cas nous avons déjà interrogé le serveur cartographique depuis le service WFS afin de disposer des éléments de représentation des départements, pourrait nous suffire dans la vue générale de l’information de décompte des établissements par département. Le décompte pourrait très bien être affiché par rapport à la géométrie des départements déjà connue.

Script de proxy à la base de données PostgreSQL

Le script distant agit comme un proxy, il permet de masquer les informations de connexion à la base de données et isoler les traitements sur la base. Pour l’application le script sera développé en langage PHP, mais évidemment n’importe quel autre langage au travers d’un serveur Web peut faire l’affaire.

Afin de simplifier le prototype j’utiliserai un script par requête spatiale. Sur un projet d’envergure on peut évidemment préférer mettre en place un contrôleur frontal.

Le premier script va effectuer une requête qui permettra de représenter le nombre d’établissements à proximité de cultures pour chaque département lorsque le niveau d’agrandissement de la carte laisse apercevoir une étendue géographique relativement large, de l’ordre de dimension d’un quart de la France métropolitaine.

Le rendu sera le suivant

Représentation par cercles proportionnels

L’analyse par symboles proportionnels sert à représenter un indicateur quantitatif en valeurs absolues (nombres, quantités, surfaces…). Ici je ne peux pas afficher en temps réel les 60 000 établissements sur la carte de France, aussi j’utiliserai des cercles proportionnels pour donner une première indication d’un nombre d’établissements adjacents aux cultures pour la distance renseignée dans le formulaire de paramétrage de la carte, pour chaque département. Cette première représentation est minimaliste et peu utile – c’est plutôt le savoir faire qui m’intéresse ici-, elle ne prend pas en compte les autres paramètres afin de pouvoir représenter un minimum d’informations en temps réel.

Chaque valeur de l’indicateur est représentée par un symbole dont la surface est proportionnelle à la valeur représentée, notre indicateur sera le décompte d’établissements adjacents.

Ces symboles proportionnels sont ajoutés sur la carte par ordre décroissant de sorte que les symboles les plus petits soient au-dessus des plus gros.

Interroger POSTGIS en PHP et retourner un résultat JSON / GeoJSON

Le script est enregistré à la racine du répertoire web du serveur Apache, il se nomme action1.php

<?php
if (!isset($_GET['callback'])) {
  header('status: 400 Bad Request', true, 400);
  exit;
}

$conn = new PDO('pgsql:host=localhost;dbname=Pesticides', 'user', 'password');
    
$distance = isset($_GET['distance']) ? (int) $_GET['distance'] : 250;

$sth = $conn->prepare('
    select p.intersections, p.code_dept, d.nom_dept, 
    ST_AsGeoJSON(ST_Centroid(d.geom)) as centroid
    from proximites p
    inner join "Departement" d
    ON p.code_dept = d.code_dept 
    where proximite = :proximite  
    order by p.intersections desc'
);
$sth->execute(array(':proximite' => $distance));
$result = $sth->fetchAll(PDO::FETCH_OBJ);

header('content-type: application/javascript; charset=utf-8');
header("access-control-allow-origin: *");
echo filter_var($_GET['callback'], FILTER_SANITIZE_ENCODED), '(', json_encode($result), ');';

Que fait ce script ?

Premièrement il vérifie la présence dans les paramètres d’url d’un paramètre nommé callback. Ce paramètre est simplement le nom de fonction qui encapsule l’objet JSON retourné (format JSONP). J’ai déjà expliqué ce principe, vous pouvez vous référer aux parties 7 ou 8 de ce document.

Le script ouvre une connexion locale sur le serveur PostgreSQL hébergé sur la même machine, puis il exécute une requête préparée sur la connexion ouverte avec le paramètre de distance fourni afin de récupérer les informations qui serviront à la représentation. Les informations de requêtes sont récupérées sous la forme d’un objet PHP grâce au paramètre PDO::FETCH_OBJ qui sera intégralement converti au format JSON par la fonction PHP json_encode().

Finalement le script renvoi les entêtes de contenu javascript et celles d’autorisation d’accès depuis un hôte distant, puis affiche le contenu de retour JSON.

Si l’on s’attarde sur la requête spatiale, vous pouvez constater que les champs code département, nom du département, et décompte sont retournés tels quels en valeurs de chaînes de caractères ou entiers.

Pour la représentation je récupère une géométrie sous forme de POINT, le centre du polygone qui défini un département.

ST_AsGeoJSON(ST_Centroid(d.geom)) as centroid

ST_Centroid – Retourne le centre géométrique d’une entité géométrie, dans le cas du département il retourne un point aux coordonnées du centre de masse de la géométrie (barycentre 2D du polygone qui matérialise les contours du département)

Le point central est transformé en notation GeoJSON pour la commodité du format d’échange.

ST_AsGeoJSON — Retourne une géométrie en tant qu’élément GeoJSON. Notez que cette fonction peut prendre en argument un nombre maxdecimaldigits qui permet de préciser la précision décimale souhaitée, une précision plus faible à notamment pour intérêt d’alléger le volume de données transférées.

Afin de respecter l’ordre de représentation des cercles, la requête effectue un tri sur le décompte, les plus petits cercles seront les derniers représentés dans l’ordre du résultat de recherche.

Exploiter le résultat de requête spatiale POSTGIS depuis un appel Ajax (côté client Javascript)

Le script PHP sert de proxy, il effectue une requête sur la base de données et retourne un résultat au format JSONP où les données spatiales sont embarquées avec les autres éléments de données mais au format GeoJSON.

Côté client Javascript il suffit donc d’appeler le proxy et traiter le résultat lorsque l’interaction avec l’utilisateur requiert de rafraîchir les informations. Nous avons deux possibilités d’interactions : soit l’utilisateur agit sur les contrôles de navigation de la carte afin de changer le niveau d’agrandissement ou de se déplacer vers une autre étendue géographique, soit il agit sur le formulaire de paramétrage pour modifier les conditions de rendu (dans ce cas uniquement la distance de proximité).

Chaque événement va provoquer le rafraîchissement de la vue en cours, à cette fin on va utiliser une fonction qui permettra de réinitialiser les couches dynamiques de la carte.

var dynamicFeatures = new L.featureGroup([]).addTo(map);

Sous Leaflet nous pouvons déclarer chaque couche représentée dynamiquement sous un groupe qui permet de manipuler toutes les couches au travers d’un seul élément. Le groupe est initialisé dans le contexte global et ajouté à la carte.

function refreshView() {
    dynamicFeatures.clearLayers(); 
...
    if (map.getZoom() < 10)
        getDecompteDepartemental();
...
}

La fonction qui rafraîchit la vue en cours commence par supprimer toutes les couches du groupe. Ensuite en fonction du niveau d’agrandissement de la carte on va rechercher les informations souhaitées. L’affichage de cercles proportionnels au niveau départemental est déclenché sous un niveau de zoom à 10, au dessus de ce niveau on choisira un autre mode de représentation.

$("#settings-form").on("change submit", function(event) {
  event.preventDefault();
  refreshView();
});

JQuery permet de créer un écouteur d’événement sur le formulaire de paramètres d’identifiant #settings-form. Les événements change et submit, déclencheront le rafraîchissement des couches dynamiques lors d’un changement sur l’un des champs du formulaire ou lors de sa soumission.

map.on('load zoomend moveend', function() {
    refreshView();
});

L’API Leaflet permet de définir directement un écouteur d’événements sur la carte. Ici on déclenche le rafraîchissement de la vue lors du chargement, lors du changement de niveau d’agrandissement, et lors d’une opération de glisser déposé.

Ces premiers éléments en place, on peut implémenter la fonction qui va représenter les cercles proportionnels.

function getDecompteDepartemental() {

    var distance = $("#distance").val();

    var query = host + '/action1.php?callback=getJson&distance=' 
            + encodeURIComponent(distance);

    $.ajax({
        url: query,
        dataType: 'jsonp',
        jsonpCallback: 'getJson',
        success: function (response) {
...
        }
    });
}

Le mécanisme est toujours le même (voir les documents antérieurs pour adaptation à la librairie OpenLayers 3), premièrement je construis l’url vers le script de proxy, cette url accepte deux paramètres, le nom de la fonction de rappel exploitée par le format JSONP et la distance récupérée en JQuery depuis le champ de sélection du formulaire qui possède l’identifiant #distance.

Un appel asynchrone permet d’obtenir la réponse du script de proxy, et en cas de succès on peut procéder à l’exploitation des résultats dans le format de réponse pour représenter les couches de données sur la carte.

En cas de succès, nous avons besoin de la valeur minimale (ou maximale) parmi les valeurs de décompte, ce minimum est exploité pour le calcul du rayon de chaque cercle.

...
            var values = [];
            for (var key in response)
                values.push(parseInt(response[key].intersections));

            var min = Math.min.apply(null, values);
...

Pour représenter chaque cercle il faut parcourir le résultat de requête, et pour chaque ligne de résultat créer un cercle proportionnel au décompte ayant pour centre le point (la géométrie) retournée au format GeoJSON qui correspond au centre de masse du polygone départemental. Ce cercle est une couche vectorielle qui est ajoutée au groupe d’éléments dynamiques pour être représentée sur la carte.

...
            for (var key in response) {
                var obj = response[key];
                var geometry = JSON.parse(obj.centroid);
                var circle = L.circleMarker([geometry.coordinates[1], geometry.coordinates[0]], {
                    color: 'white',
                    fillColor: 'Orange',
                    fillOpacity: 0.25,
                    opacity: 1
                });
                circle.setRadius(getRadius(parseInt(obj.intersections), min));
                circle.bindPopup(obj.nom_dept + " : " 
                        + obj.intersections + " établissement(s) à proximité" );
                
                dynamicFeatures.addLayer(circle);
            }
...

Pour chaque élément de réponse, on assigne l’élément à un objet nommé obj. Cet objet contient les propriétés telles que définies dans notre requête SQL, c’est à dire : intersections, code_dept, nom_dept, centroid. Les intersections sont un nombre, le décompte, le code département et le nom de département sont des chaînes de caractères, et le centroid un point au format GeoJSON, le format GeoJSON est parcouru pour être transformé en objet javascript à l’aide de la fonction Javascript JSON.parse(). Depuis ce nouvel objet on peut récupérer les coordonnées de latitude et longitude du point pour la représentation sous forme de cercle.

Le rayon du cercle est défini par rapport au décompte et au décompte minimal puis un popup est rattaché au cercle, il s’affichera au clic sur le cercle. Le contenu du popup est au format html, il contient des informations relatives au départements récupérées depuis le résultat de requête.

Dernier point crucial : le cercle est ajouté au groupe que l’on a créé en amont et qui est déjà inclus dans la carte, ce qui provoque l’affichage de l’élément sur la carte.

Concernant le calcul du rayon pour la représentation du cercle je vous communique les fonctions sans entrer dans les détails. Je ne suis pas entièrement satisfait du résultat lors de changements d’échelles de la carte et je vous invite à proposer une solution améliorée.

...
// calcule le rayon du cercle proportionnel
// Retourne le rayon du plus petit cercle

function getMinRadius() {
    var bounds = map.getPixelBounds();
    return Math.floor((bounds.max.x - bounds.min.x) * 0.005);
//    var size = map.getSize();
//    return Math.floor(size.x * 0.005);
}

// Retourne le rayon du cercle
// n est la valeur
// min est la valeur minimale

function getRadius(n, min) {
    return Math.floor(getMinRadius() * Math.sqrt(n / min));
}
...

Sur le même principe nous pouvons récupérer et représenter des géométries plus complexes, c’est exactement ce dont on a besoin à un niveau d’agrandissement élevé pour représenter un établissement.

Un établissement est représenté par un marqueur de position, une icône, la distance de proximité à l’établissement est matérialisée par un cercle opaque autour de ce point, et l’intersection entre les cultures et ce disque est mise en évidence par une couche opaque superposée.

Le script de proxy, nommé action2.php est très similaire au script précédent, seul le traitement des paramètres d’url et la requête à la base changent.

La requête SQL est un cas très interessant, riche en enseignements, et c’est le coeur du projet aussi nous allons la passer en revue. Je vais ignorer le traitement en amont des paramètres pour être plus succint.

...
$bbox = explode(',', $_GET['bbox']);

$sql = '
WITH resultats AS (
SELECT 
e.numero_uai, e.appellatio, e.adresse_ua, e.code_post, e.localite_a, e.nature_uai,
ST_AsGeoJSON(ST_Transform((e.geom), 4326), 6) AS etablissement,
floor(sum(ST_Area(CAST(cul_geom As geography)))) * POWER(10, -4) as surface_totale_cultures,
floor(sum(ST_Area(CAST(ST_Intersection(ST_Buffer(e.geom::geography, :distance)::geometry, cul_geom) As geography)))) * POWER(10, -4) as surface_culture_proximite,
ST_AsGeoJSON(ST_Transform((ST_Union(cul_geom)), 4326), 6) AS cultures,
ST_AsGeoJSON(ST_Transform((ST_Union(ST_Intersection(ST_Buffer(e.geom::geography, :distance)::geometry, cul_geom))), 4326), 6) AS surfaces_proximite,
count(i.num_ilot) as nombre_parcelles,
count(distinct c.cult_maj) as types_culture 
FROM intersect_etabs i 
INNER JOIN "Etablissement" e 
ON i.numero_uai = e.numero_uai 
INNER JOIN master_cultures c
ON i.num_ilot = c.num_ilot  
WHERE cul_geom && ST_MakeEnvelope(:xmin, :ymin, :xmax, :ymax, 4326)
AND distance < :distance 
AND c.cult_maj not in (18, 19)  -- Elimine les prairies
AND e.nature_uai < 800          -- Elimine les etablissements administratifs
AND e.etat_etabl = 1            -- Elimine les etablissements fermes
' . $whereClause . ' 
group by e.numero_uai  
order by e.numero_uai
)
SELECT *
FROM resultats
WHERE surface_culture_proximite > :surface
' . $aggregateWhereClause . ' 
ORDER BY types_culture DESC, surface_culture_proximite DESC, nombre_parcelles DESC
';

$sth = $conn->prepare($sql);
$sth->execute(array(
    ':distance' => $distance,
    ':surface' => ($surface * pow(10, -4)),
    ':xmin' => $bbox[0],
    ':ymin' => $bbox[1],
    ':xmax' => $bbox[2],
    ':ymax' => $bbox[3],
));
...

Cette requête est assez monolithique à première vue, nous pouvons la découper pour appréhender le fonctionnement.

...
WITH resultats AS (
SELECT 
...
ST_AsGeoJSON(ST_Transform((ST_Union(ST_Intersection(ST_Buffer(e.geom::geography, :distance)::geometry, cul_geom))), 4326), 6) AS surfaces_proximite,
...
group by e.numero_uai
)
SELECT *
FROM resultats
WHERE surface_culture_proximite > :surface
...

L’application permet de sélectionner une surface dans le formulaire de paramètres, l’un des objectifs est donc de déterminer la surface totale impactée dans le rayon de proximité autour de l’établissement toutes parcelles confondues pour appliquer un filtre par rapport au paramètre.

Pour calculer la surface il est nécessaire d’agréger les données de parcelles pour chaque établissement. L’opération d’agrégat est réalisée par l’instruction SQL GROUP BY, le résultat d’une requête de type SELECT va retourner les résultats agrégés, une autre requête est alors exécutée sur le résultat de la précédente pour appliquer les paramètres de filtre.

L’imbrication de requêtes est réalisée par les instructions WITH, AS de PostgreSQL, les clauses de conditions stipulent de retourner les résultats agrégés pour lesquels la surface est supérieure à la surface paramétrée (c’est également au niveau de ces clauses SQL que l’on insère si besoin la restriction d’adjacence à plusieurs cultures)

Comment aggréger des objets géométriques ?

Dans notre cas on définit en premier une zone tampon autour du point qui matérialise l’établissement avec la fonction ST_Buffer, que nous avons déjà employé précédemment.

ST_Buffer(e.geom::geography, :distance)::geometry

ST_Intersection — Retourne une géométrie qui représente la portion commune de deux géométries, ici elle est utilisée sur la géométrie de la zone tampon et celle de la parcelle pour ne conserver que la surface en commun.

ST_Union — Retourne une géométrie qui représente le jeu de points d’union ensembliste des géométries, elle va réunifier toutes les intersections de la zone tampon avec les cultures à proximité de l’établissement.

ST_Transform – Retourne une nouvelle géométrie avec ses coordonnées transformées pour le SRID référencé par le paramètre entier, dans notre cas le SRID 4326 permet de travailler selon nos attentes depuis la librairie Javascript.

Une fois ces informations comprises, le reste de la requête est presque trivial. Expliquons les deux derniers points de difficulté.

SELECT ... floor(sum(ST_Area(CAST(cul_geom As geography)))) * POWER(10, -4) as surface_totale_cultures,

L’information retournée comporte une indication de la surface totale des parcelles à partir du moment où la parcelle répond aux paramètres de recherche. La surface est calculée depuis un objet géographie.

L’opération de conversion de type permet de transformer la géométrie de la parcelle en géographie, la fonction ST_Area va retourner l’aire de l’objet, tandis que la fonction SUM va additionner les aires.

Au final la somme des aires est arrondie à l’entier le plus bas et multipliée par la puissance de 10 à l’exposant -4 afin de retourner une quantité exprimée en hectares.

La fonction ST_MakeEnvelope (voir plus haut) restreint la requête aux résultats sur l’étendue géographique passée en paramètres depuis le client Javascript.

    var query = host + '/action2.php?callback=getJson&' + $("#settings-form").serialize()
        + '&bbox=' + bounds.toBBoxString();

Les paramètres sont construits depuis la fonction d’appel, et le résultat de requête est traité à nouveau en cas de succès.

L’Education Nationale a mis à disposition un répertoire des établissements scolaires, l’application permet de renseigner le numéro d’établissement dans le formulaire de recherche.

Pour clôre ce chapitre, voici côté client la fonction de représentation des établissements, suivie d’une copie d’écran du rendu sur la carte.

var premierDegre = L.MakiMarkers.icon({icon: "school", color: "#b0b", size: "m"});
var secondDegre = L.MakiMarkers.icon({icon: "school", color: "#E2492F", size: "m"});

function getEtablissements(bounds) {

    var distance = $("#distance").val();

    var query = host + '/action2.php?callback=getJson&' + $("#settings-form").serialize()
        + '&bbox=' + bounds.toBBoxString();

    $.ajax({
        url: query,
        dataType: 'jsonp',
        jsonpCallback: 'getJson',
        success: function (response) {

            for (var key in response) {
                var obj = response[key];

                // Représente l'établissement par un marqueur de position
                
                var etablissement = JSON.parse(obj.etablissement);
                var coordinates = etablissement.coordinates[0];
                var icon = (obj.nature_uai < 300) ? premierDegre : secondDegre;
                var link = "<a href='http://www.education.gouv.fr/bce/index.php?simple_public=" 
                        + obj.numero_uai + "'>" + obj.appellatio + "</a>"
                
                var marker = L.marker([coordinates[1], coordinates[0]], {icon: icon})
                    .bindPopup(
                        link + "<br /> " 
                        + (obj.adresse_ua ? obj.adresse_ua + ', ' : '') 
                        + obj.code_post + ' ' + obj.localite_a + "<br /> " 
                        + 'Surface à proximité : ' + obj.surface_culture_proximite  + " ha<br /> " 
                        + 'Surface totale parcelles : ' + obj.surface_totale_cultures + " ha<br />"
                        + obj.nombre_parcelles + ' parcelle(s)' + "<br />"
                        + obj.types_culture + ' type(s) de cultures' + "<br />"
                        );

                // La zone tampon autour de l'établissement est représentée par un cercle

                var buffer = L.circle([coordinates[1], coordinates[0]], distance, {
                    color: 'red',
                    fillColor: '#f03',
                    fillOpacity: 0.25,
                    opacity: 0.65,
                    weight: 1
                });

                var parcelles = L.geoJson(JSON.parse(obj.cultures), {
                    style: {
                        "color": "#ff7800",
                        "weight": 1,
                        "opacity": 0.65
                    }
                });
                
                var intersections = L.geoJson(JSON.parse(obj.surfaces_proximite), {
                    style: {
                        "color": "red",
                        "weight": 1,
                        "opacity": 0.95
                    }
                });
                
                dynamicFeatures.addLayer(parcelles);
                dynamicFeatures.addLayer(intersections);
                dynamicFeatures.addLayer(buffer);
                dynamicFeatures.addLayer(marker);   
                
                if (map.getZoom() > 15)
                    getJsonLayerParcelles(bounds);
            }
        }
    });
}

Fonctionnalités annexes spécifiques Leaflet

Pour perfectionner un peu la carte, j’ai ajouté un ensemble de fonctionnalités.

Lorsque l’utilisateur arrive sur la carte, j’ai trouvé utile de positionner la carte sur sa localisation.

Cela peut se faire nativement sous Leaflet.

map.locate({setView: true, maxZoom: 15});

La barre latérale est par défaut repliée, ce qui n’est pas intuitif pour l’utilisateur mais si l’on déplie la barre par défaut il n’est pas intuitif de la replier non plus.

En attendant un boutton de fermeture sur cette extension, j’ouvre la barre latérale au chargement et je la referme après quelques secondes. Les opérations sur la barre peuvent s’effectuer depuis les identifiants html des onglets.

setTimeout(function () {
    sidebar.open('home');
}, 1000);
setTimeout(function () {
    sidebar.close('home');
}, 4000);

Finalement, afin de donner à l’utilisateur la possibilité de passer en vue aérienne je fais appel aux serveurs WMS de l’IGN.

Conclusion

Au terme de cette série de documents sur la réalisation d’une application cartographique j’ai détaillé comment réaliser une carte dynamique depuis les données et leur mise à disposition, jusqu’à la représentation dynamique des données.

Techniquement je n’ai pas abordé par manque de temps la représentation d’une carte par chloroplèthe, ce sera l’objet d’un futur document, cette fois ci sous OpenLayers 3.

J’espère que cette application trouvera une utilité publique ou qu’à défaut ce travail motivera ce type de représentations.

Requête à un serveur cartographique sous client web OpenLayers3

Etude du code prototype projet OpenLayers 3

Nous avons généré automatiquement dans QGIS un prototype qui permet de visualiser une carte au travers d’un client web, explorons le code…

Le meilleur moyen d’aborder OpenLayers3 est de passer par les exemples et d’explorer la librairie. (Téléchargez la dernière version d’OpenLayers, vous trouverez les codes sources des exemples et la librairie en plus de la version des codes sources au complet)

Concept

La composante de base d’OpenLayers3 est la carte (ol.Map). Cette carte est rendue dans un conteneur cible, en pratique dans un élément div de la page Web qui contient la carte.

Les propriétés de la carte peuvent être configurées au moment de sa construction, ou ultérieurement par des méthodes spécifiques.

Le conteneur cible qui a été généré pour nous est un élément div d’identifiant map (il comporte des sous-divisions pour permettre un affichage en fenêtre modale dont on pourrait se passer pour simplifier l’exemple…).

        <div id="map">
            <div id="popup" class="ol-popup">
                <a href="#" id="popup-closer" class="ol-popup-closer"></a>
                <div id="popup-content"></div>
            </div>  
        </div>

Ce conteneur possède un style d’affichage paramétré en entête du code html, le style est paramétré ici pour un affichage à 100% de la taille de la fenêtre du navigateur.

        <style>
            html, body, #map {
                height: 100%;
                width: 100%;
                overflow: hidden;
            }
         ...
        </style>

Notez que la librairie OpenLayers et ses extensions possèdent leurs propres fichiers de style qui sont chargés en ressources en entête de fichier html.

        <link rel="stylesheet" href="./resources/ol.css">
        <link rel="stylesheet" href="./resources/boundless.css">

ol.Map ne gère pas lui même le niveau de zoom et la projection de la carte, ce sont des propriétés et méthodes d’une instance de vue, ol.View, qui est initialisée dans la carte.

Pour obtenir les données à représenter pour une couche, OpenLayers utilise des sous-classes ol.source.Source. Ces sources permettent d’alimenter la carte au travers de serveurs de tuiles, de serveurs WMS ou WMTS, de données vectorielles dans des formats comme GeoJSON ou KML. Voyons comment sont configurées les sources par défaut dans notre prototype.

Chargement des ressources javascript

Premièrement le prototype charge la librairie OpenLayers 3

<script src="./resources/ol.js"></script>

Puis une il fait appel à une extension qui implémente un composant qui permet de permuter la visibilité des couches…

<script src="./resources/boundless.js"></script>

Le script charge ensuite les données de chaque couche, ces données ont été transformées au format Geojson par l’extension installée sous QGIS, ce sont ces données qui sont les sources de notre prototype par défaut.

Si l’on considère la couche qui comporte les informations départementales, le script va faire appel à deux ressources, la première comporte les données de la couche au format Géojson

<script src="layers/Departement.js"></script>

qui a pour contenu

var geojson_Departement = {
    "type":"FeatureCollection"
    },
    "crs":{
    "type":"name",
            "properties":{
            "name":"urn:ogc:def:crs:EPSG::3857"}
    }, "features":[
    {"type":"Feature", "properties":{}, "geometry":{"type":"MultiPolygon", "coordinates":
            [[[[649129.155, 5770492.803], [648115.547, 5769160.517],
                    [648898.518, 5767096.316], [648562.967, 5764862.309],
                    [647339.990, 5760692.349], [646490.730, 5758044.713], 
                            ...

GeoJson est un format qui permet d’encoder des structures de données géographiques, nous ne nous attarderons pas sur le format, notre objectif est de remplacer ces données statiques par des données fournies dynamiquement par notre serveur cartographique.

La seconde ressource est un script Javascript qui configure le style de représentation de la couche.

<script src="styles/Departement_style.js"></script>

qui a pour contenu (relativement obscure à cause du code environnant)

var styleCache_Departement = {}
var style_Departement = function (feature, resolution) {
    var value = ""
    var style = [new ol.style.Style({
            stroke: new ol.style.Stroke({color: "rgba(0,255,0,1.0)", lineDash: null, width: 1})
        })
    ];
    var labelText = "";
    var key = value + "_" + labelText

    if (!styleCache_Departement[key]) {
        var text = new ol.style.Text({
            font: '8.25px Calibri,sans-serif',
            text: labelText,
            fill: new ol.style.Fill({
                color: "rgba(0, 0, 0, 255)"
            }),
        });
        styleCache_Departement[key] = new ol.style.Style({"text": text});
    }
    var allStyles = [styleCache_Departement[key]];
    allStyles.push.apply(allStyles, style);
    return allStyles;
};

Finalement le programme charge une dernière ressource qui permet d’initialiser chaque couche depuis les ressources chargées précédemment.

<script src="./layers/layers.js" type="text/javascript"></script>

C’est dans ce script que se passe le plus gros travail d’initialisation.

var baseLayer = new ol.layer.Tile({title: 'OSM', source: new ol.source.OSM()});
...
var lyr_Departement = new ol.layer.Vector({
                source: new ol.source.GeoJSON({object: geojson_Departement}), 
                style: style_Departement,
                title: "Departement"
            });
...
lyr_Departement.setVisible(true);
...
var layersList = [baseLayer,lyr_cultures,lyr_Departement,lyr_DepartementWFS,lyr_etabstamponselec,lyr_etabssel,lyr_etabparcelintersm];

Initialisation des couches

Dans le script layers.js, le script va construire la liste des couches à afficher. Chaque couche est crée l’une à la suite de l’autre avec une source de données, un style et d’autres paramètres. Puis la visibilité de chaque couche est positionnée à la valeur vraie (c’est toutefois la valeur par défaut lors de l’initialisation). La liste des couches est le tableau résultant.

var baseLayer = new ol.layer.Tile({title: 'OSM', source: new ol.source.OSM()});

Ici on construit un objet ol.layer.Tile qui va permettre d’interroger un fournisseur de tuiles (images pré-rendues). L’argument title n’est pas dans les options de l’objet, mais comme cet objet implémente la classe abstraite ol.Object l’argument permet de renseigner un attribut qui possédera des accesseurs (qui pourront être exploités par le composant qui permet de permuter la visibilité des couches)

La source de données ol.source.OSM est pré-configurée dans la librairie et permet d’interroger les serveurs de tuiles d’OpenStreetMap.

var lyr_Departement = new ol.layer.Vector({
                source: new ol.source.GeoJSON({object: geojson_Departement}), 
                style: style_Departement,
                title: "Departement"
            });

Ces instructions permettent de créer une couche de données vecteur, rendue côté client. La source de la couche est au format GeoJSON, elle accepte un objet GeoJSON, geojson_Departement, qui a été construit lors de l’export du projet QGIS et chargé en amont dans les ressources Javascript.

Le style de rendu de la couche vectorielle ol.layer.Vector est également fournit en option lors de la construction, il permet de passer l’objet ol.style.Style également chargé en amont dans les ressources javascript.

lyr_Departement.setVisible(true);

L’objet de couche est déjà défini, le script utilise explicitement un accesseur pour conditionner la visibilité.

var layersList = [baseLayer, ..., lyr_Departement, ...];

La liste des couches est finalement renseignée dans une variable layersList a contexte global.

Initialisation et rendu de la carte

Le script vient d’obtenir les données à représenter pour chaque couche il reste à initialiser l’objet carte mais avant cela il faut encore expliquer un morceau de code généré pour nous.

            var container = document.getElementById('popup');
...
            var overlayPopup = new ol.Overlay({
                element: container
            });

Sous OpenLayers3 il existe une classe de widgets qui s’affiche en superposition, ol.Overlay.

Le code créé une instance de widget qui sera superposée à l’élément du DOM récupéré dans la variable container (qui possède l’identifiant popup)

L’extension disponible sous QGIS génère beaucoup de code superflus pour un usage basique. La gestion de popups et overlays ne m’est d’aucune utilité et sera supprimée du prototype ainsi que toute la gestion d’événements à la souris.
            var map = new ol.Map({
                controls: ol.control.defaults().extend([
                    new ol.control.ScaleLine({}), new Boundless.LayersControl({groups: {default: {title: "Layers"}}})
                ]),
                target: document.getElementById('map'),
                renderer: 'canvas',
                overlays: [overlayPopup],
                layers: layersList,
                view: new ol.View2D({
                    maxZoom: 28, minZoom: 1
                })
            });

La construction de l’objet carte ol.Map est le cœur de la librairie. Elle requiert de renseigner une vue, ajouter une ou plusieurs couches, et préciser un élément conteneur cible.

Les options générées dans l’exemple sont les suivantes :

  • controls : Collection de contrôles. Un contrôle ol.control.Control est un widget représenté en position fixe à l’écran. Un contrôle peut être seulement informatif ou bien peut permettre d’interagir avec l’utilisateur. Une collection de contrôles par défaut est incluse aux cartes, ol.control.defaults (Elle comporte les widgets de Zoom, Rotation et Attributions). Ici le code ajoute deux contrôles à la collection par défaut, en premier l’échelle et en second le widget qui permet de sélectionner les couches à afficher ou masquer. Le rendu des contrôles est paramétrable par application de styles css
  • target : Elément du DOM où la carte sera incluse
  • renderer : Moteur de rendu, seul le canvas a les capacités de rendu vectoriel
  • overlays : Collection de widgets en superposition (contrairement aux contrôles la position des widget n’est pas fixe, elle dépend de la carte)
  • layers : Collection de couches
  • view : Vue associée à la carte.

La vue est un objet qui représente une vue 2D simple de la carte, elle permet de spécifier le centre, la résolution et la rotation de la carte. La façon la plus simple de définir une vue est de définir un point central et un niveau de zoom.

Dans notre exemple la vue est initialisée uniquement avec un niveau de zoom minimal et maximal.

Un objet ol.View possède une projection. La projection détermine le système de coordonnées du centre et ses unités déterminent les unités de résolution (projection d’unités par pixels). La projection par défaut est sphérique Mercator (EPSG: 3857), système de projection utilisé entre autre par OSM et Google Maps, et comme vous pouvez constater elle n’est pas renseignée dans notre cas et correspond bien à la projection utilisée dans les fichiers de données Geojson.

Et ensuite ?

map.getView().fitExtent([-225069.813761, 5267288.926421, 186886.432518, 5746435.175340], map.getSize());

La vue associée à l’objet carte est récupérée. La méthode fitExtent() est au stade expérimental (vous devez décochez la case Stable Only dans l’API pour faire apparaître la documentation !), elle fait s’adapter la vue à l’étendue géographique et à la taille données. La taille est en pixels, elle est déterminée ici par la taille de la carte. L’étendue est un tableau de nombres qui représentent une étendue géographique, il s’agit de [minx, miny, maxx, maxy]. Ces valeurs sont les valeurs récupérées sous QGIS lors de l’export du projet, il est bien entendu possible de les modifier pour s’adapter à une autre étendue par défaut. (Il suffit sous QGIS de copier-coller les valeurs pour la vue qui convient)

Comme les données ne concernent que la Gironde, j’ai basculé le SCR du projet à EPSG:3857 sous QGIS et zoomé la sélection sur la Gironde à partir des attributs de la couche des départements. L’emprise est disponible dans la barre de statut de QGIS en bas à gauche après permutation de l’affichage de position du curseur.

Ce qui donne donc

map.getView().fitExtent([-174152,5464320,85720,5738798], map.getSize());

Le reste du code n’a pas d’incidence sur le rendu de carte, il permet de gérer l’interaction souris avec l’utilisateur. (Et à vrai dire ne semble pas fonctionner)

Nous avons donc expliqué le code et vu comment générer une carte à partir de sources de données vectorielles sous forme de fichiers. Comment adapter le code pour interroger un serveur de données ?

Adaptation du code pour une source de données par serveur cartographique

Nous avons a disposition pour le projet deux services qui tournent sous un serveur cartographique QGIS Server.

Les capacités WFS et WMS des services peuvent être interrogées par les url suivantes

http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetCapabilities
http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

C’est ce dernier service que l’on va exploiter en premier.

Appel à un service WMS QGIS Server sous OpenLayers3

Le nom de la couche des départements nous intéresse particulièrement, on peut le récupérer depuis le fichier xml.

<Layer queryable="1">
<Name>Departement</Name>
<Title>Departement</Title>
<CRS>EPSG:4326</CRS>
<EX_GeographicBoundingBox>
<westBoundLongitude>-5.21251</westBoundLongitude>
<eastBoundLongitude>9.63332</eastBoundLongitude>
<southBoundLatitude>41.3141</southBoundLatitude>
<northBoundLatitude>51.1376</northBoundLatitude>
</EX_GeographicBoundingBox>
<BoundingBox CRS="EPSG:4326" maxx="51.1376" minx="41.3141" maxy="9.63332" miny="-5.21251"/>
<Style>
<Name>default</Name>
<Title>default</Title>
<LegendURL>
<Format>image/png</Format>
<OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetLegendGraphic&LAYER=Departement&FORMAT=image/png&STYLE=default&SLD_VERSION=1.1.0"/>
</LegendURL>
</Style>
</Layer>

Le parcours des sources de données ol.source dans l’API OpenLayers3 nous renseigne sur les possibilités de la librairie.

Il apparaît possible d’utiliser un service WMS pour fournir les images pour une couche donnée, le type est ol.source.ImageWMS.

Il suffit de remplacer la source locale pour la couche des départements.

//var lyr_Departement = new ol.layer.Vector({
//                source: new ol.source.GeoJSON({object: geojson_Departement}), 
//                style: style_Departement,
//                title: "Departement"
//            });

var lyr_Departement = new ol.layer.Image({
                source: new ol.source.ImageWMS({
                               url: "http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi",
                               params: {
                                   'LAYERS': 'Departement'
                               },
                            }),
                style: style_Departement,
                title: "Departement"
            });

La source ol.source.ImageWMS requiert l’url du service et un tableau des paramètres de la requête WMS en options, dont seul le paramètre layers est obligatoire. (Les autres paramètres sont renseignés dynamiquement par la librairie)

Il est possible de faire la même manipulation pour les autres couches en modifiant le nom de la couche dans l’option layers.

Les requêtes vers le serveur WMS peuvent être contrôlées sous les outils de développement de votre navigateur Web préféré afin d’étudier la requête générée et le temps d’exécution. A titre d’exemple une requête effectuée pour récupérer les parcelles de culture.

http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&FORMAT=image%2Fpng&TRANSPARENT=true&LAYERS=cultures_33&CRS=EPSG%3A3857&STYLES=&WIDTH=2880&HEIGHT=1529&BBOX=-33130.75605296711%2C5560982.287734281%2C-5613.425870303659%2C5575591.314765286

Notez également que depuis ce type d’appel les paramètres d’affichage en résolution minimale et maximale sont respectés.

var lyr_cultures = new ol.layer.Image({
                source: new ol.source.ImageWMS({
                               url: "http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi",
                               params: {
                                   'LAYERS': 'cultures_33'
                               },
                            }),
                minResolution:-4.65661009752e-10,
                maxResolution:55785.0,
                style: style_cultures,
                title: "cultures_33"
            });

Dans l’exemple précédent les cultures ne s’affichent que dans la résolution souhaitée par contre les requêtes sont effectuées sur le serveur quelle que soit la résolution dans la vue représentée. Hors limites de représentation elles retournent une image vide (par défaut un png transparent).

Afin de limiter le nombre de requêtes au serveur il peut être intéressant de demander plusieurs couches dans une même requête, par exemple on peut regrouper les trois couches qui servent à représenter les établissements et leur proximité aux cultures.

var lyr_etabs = new ol.layer.Image({
        source: new ol.source.ImageWMS({
            url: "http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi",
            params: {
                'LAYERS': 'etabs_tampon_selec,etabs_sel,etab_parcel_inters_50m'
            },
        }),
        title: "Etablissements"
    });

En passant par le service WMS le serveur cartographique retourne une image, de ce fait on perd les possibilités offertes par la couche vectorielle Geojson. L’avantage d’une image est souvent sa légèreté sur le réseau par rapport à un fichier de données, mais le format vectoriel nous permet d’interagir en javascript avec les objets, on peut par exemple sélectionner un objet et modifier son style, supprimer ou ajouter un objet.

Interagir avec les objets peut permettre de créer des cartes plus interactives, en créant des composants de formulaire en widget par exemple on peut imaginer de représenter dynamiquement une carte choroplèthe ou des symboles proportionnels en fonction d’un paramètre sélectionné et d’une requête asynchrone qui va chercher les valeurs d’indicateurs à représenter.

Appel à un service WFS QGIS Server sous OpenLayers3

La librairie OpenLayers 3 a la capacité de charger des couches vectorielles depuis une méthode ol.layer.Vector() comme nous l’avons vu précédemment. A plus bas niveau elle permet d’utiliser la méthode ol.source.ServerVector(), toutefois encore en état expérimental, qui va permettre d’alimenter une source dans un format supporté depuis un serveur distant.

Il nous suffit en théorie d’utiliser le format ol.format.GeoJSON() avec notre serveur QGIS Server depuis une requête asynchrone côté client web.

En pratique, si l’on effectue une requête sur notre serveur afin de récupérer des données JSON depuis une page html côté client (en Ajax donc), le navigateur va interdire l’opération si le serveur n’autorise pas explicitement les accès en origine croisée par l’ajout d’une entête de contrôle d’accès.

Notre problème est double. Premièrement QGIS Server possède les capacités pour servir du GeoJSON, qui est simplement du JSON avec une convention, par contre il ne propose pas le format de sortie JSONP qui encapsulerait simplement le GeoJSON dans une fonction Javascript. Il nous resterait comme option de configurer le serveur pour accepter les requêtes d’origine croisées (CORS), cependant le serveur QGIS fonctionne en mode CGI et Apache ne peut pas à ma connaissance ajouter une entête d’autorisation sur la sortie produite par QGIS Server.

La solution lorsque l’on ne peut ni modifier le format de sortie ni configurer le serveur pour ajouter des autorisations est d’utiliser un proxy. Très simplement un script sera hébergé sur le même serveur, il acceptera en paramètres un nom de fonction qui servira à l’encapsulation javascript du format JSON, et l’url qu’il sera en mesure d’interroger depuis le même domaine.

J’ai créé un simple script en langage PHP qui servira de proxy. Vous pourrez probablement trouver des solutions plus complètes…

<?php
if(!isset($_GET['callback']) || !isset($_GET['url'])) {
  header('status: 400 Bad Request', true, 400);
  exit;
}
extract(parse_url($_GET['url']));       

if (! ( in_array($host, array($_SERVER['HTTP_HOST'], $_SERVER['SERVER_ADDR']) )
      || ( basename($path) != 'qgis_mapserv.fcgi' ) )) {
  header('status: 400 Bad Request', true, 400);
  exit;
}
if (!$content = file_get_contents(filter_var($_GET['url'], FILTER_SANITIZE_URL))) {
  header('status: 400 Bad Request', true, 400);
  exit;
}
header('content-type: application/javascript; charset=utf-8');
header("access-control-allow-origin: *");
echo filter_var($_GET['callback'], FILTER_SANITIZE_ENCODED), '(', $content, ');';

Ce script nommé proxy.php est placé à la racine du répertoire web. Il doit être appelé avec deux paramètres, callback et url.

Le paramètre url doit être encodé…

http://io.gchatelier.fr/proxy.php?callback=getDepartementJson&url=...

Le script effectue quelques contrôles, premièrement il n’aboutit pas si les paramètres ne sont pas renseignés. Il vérifie ensuite que l’url à interroger est sur le même domaine que le serveur, puis vérifie que l’on cherche bien à interroger un script CGI qui correspond à QGIS Server. Le script filtre ensuite les paramètres et exécute la requête. Le contenu de la requête est retourné encapsulé dans la fonction de rappel Javascript avec les bons entêtes : premièrement l’entête d’application javascript (et non pas l’entête d’application JSON), puis l’entête de contrôle d’accès qui autorise le proxy à être interrogé depuis un autre domaine.

Une fois ce problème réglé on peut passer au code.

            var DepartementWFSURL = 'http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?' + 
                    'SERVICE=WFS&REQUEST=GetFeature&' + 
                    'VERSION=1.1.0&TYPENAME=Departement&' + 
                    'SRSNAME=EPSG:3857&outputFormat=GeoJSON';
            
            var proxy = 'http://io.gchatelier.fr/proxy.php?callback=getDepartementJson&url=';

Construit l’url d’appel à notre proxy en précisant la fonction de rappel et l’url d’origine du service.

            var deptsSource = new ol.source.ServerVector({
                format: new ol.format.GeoJSON(),
                loader: function(extent, resolution, projection) {
                    var url = DepartementWFSURL +
                        '&bbox=' + extent.join(',');

                    $.ajax({
                        url: proxy + encodeURIComponent(url),
                        dataType: 'jsonp',
                        jsonpCallback: 'getDepartementJson',
                    });
                },
                strategy: function() {
                    return [ [-5.21251302, 41.31412406, 9.63331895, 51.13763146] ];
                },
                projection: 'EPSG:3857'
            });

La source récupère les données WFS en GeoJSON en utilisant la technique JSONP.

Remarquez l’option strategy. OpenLayers 3 possède plusieurs types de stratégies ol.loadingstrategy pour le chargement des données. Il peut charger les données par tuiles en multipliant les appels au serveur, tout charger d’un coup, ou charger suivant les coordonnées de la boîte englobante (BBOX). Ici les coordonnées sont fournies sur l’étendue géographique complète [minx, miny, maxx, maxy].

L’option loader correspond à la fonction qui effectue le chargement, elle se voit passer l’étendue et la projection que l’on a paramétrées. Avec la stratégie de tuiles elle serait appelée de multiples fois. C’est que l’appel javascript asynchrone va récupérer les objets de la couche vectorielle depuis le service WFS

L’appel Ajax accepte ces paramètres :

  • url : url vers notre proxy qui interrogera le serveur WFS
  • dataType : notre proxy retourne du javascript au format JSONP, on doit préciser le type de données car les entêtes retournées sont application javascript
  • jsonpCallback : le nom de la fonction de rappel dans le fichier JSONP retournée par le proxy
            var getDepartementJson = function(response) {
                deptsSource.addFeatures(deptsSource.readFeatures(response));
            };

En cas de réussite la fonction anonyme de rappel est exécutée avec la réponse en paramètre, notre GeoJSON… la réponse est parcourue et ajoutée à la source vectorielle.

            var lyr_Departement = new ol.layer.Vector({
                source: deptsSource,
                title: 'Departement',
                style: style_Departement,
            });

Il ne reste plus qu’à exploiter la source vectorielle que nous venons de créer comme nous avons exploité les sources GeoJSON.

Nous avons fait le tour des codes sources générés sous une extension QGIS pour la librairie OpenLayers3. Vous pouvez également vous référer au même document adapté pour la librairie Leaflet.

Nous sommes désormais en mesure de construire des applications interactives pour le web depuis des services cartographiques, ce sera l’objet du prochain document de cette série qui détaillera la création d’une carte interactive sous la librairie Leaflet, le même principe pouvant s’appliquer à OpenLayers3 avec les techniques que nous venons d’expliquer.

Requête à un serveur cartographique sous client web Leaflet

Etude du code prototype projet Leaflet

Nous avons généré automatiquement dans QGIS un prototype qui permet de visualiser une carte au travers d’un client web, explorons le code…

Le meilleur moyen d’aborder Leaflet est de passer par le tutoriel de démarrage et d’explorer la librairie.

Concept

La composante de base de Leaflet est la carte (L.map). Cette carte est rendue dans un conteneur cible, en pratique dans un élément div de la page Web qui contient la carte.

Les propriétés de la carte peuvent être configurées au moment de sa construction, ou ultérieurement par des méthodes spécifiques.

Le conteneur cible qui a été généré pour nous est un élément div d’identifiant map

<div id="map"></div>

Ce conteneur possède un style d’affichage paramétré dans le fichier de style css/own_style.css, le style est paramétré ici pour un affichage à 100% de la taille de la fenêtre du navigateur.

...
	html, body, #map {
		height: 100%;
		width: 100%;
		padding: 0;
		margin: 0;
	}
...

Notez que la librairie Leaflet possède son propre fichier de style qui est chargé en ressource en entête de fichier html depuis un CDN.
Les fichiers de style des extensions Leaflet.markercluster et Leaflet.label sont également chargés, ainsi que la feuille de style personnalisée pour le projet.

        <link rel="stylesheet" href="http://cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.3/leaflet.css" />
        <link rel="stylesheet" href="css/MarkerCluster.css" />
        <link rel="stylesheet" href="css/MarkerCluster.Default.css" />
        <link rel="stylesheet" type="text/css" href="css/own_style.css">
        <link rel="stylesheet" href="css/label.css" />

L.Map gère le niveau de zoom et la projection de la carte.

Pour obtenir les données à représenter pour une couche, Leaflet utilise des classes pour chaque type de couche vecteur, ou raster. Ces couches permettent d’alimenter la carte au travers de serveurs de tuiles, serveur WMS, images, données vectorielles et format GeoJSON. Voyons comment sont configurées les sources par défaut dans notre prototype.

Chargement des ressources javascript

Premièrement le prototype charge la librairie JQuery que l’on ne présente plus, puis la librairie Autolinker qui est un utilitaire qui permet de formater proprement des urls en javascript.

        <script src="http://code.jquery.com/jquery-1.11.1.min.js"></script>
        <script src="js/Autolinker.min.js"></script>

Puis une il fait appel à la librairie Leaflet et à ses extensions, leaflet-hash, Leaflet.label et Leaflet.markercluster.

        <script src="http://cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.3/leaflet.js"></script>
        <script src="js/leaflet-hash.js"></script>
        <script src="js/label.js"></script>
        <script src="js/leaflet.markercluster.js"></script>

Que font ces extensions ?

  • leaflet-hash va ajouter dynamiquement à l’url du navigateur une ancre qui permet aux utilisateurs de repérer une vue spécifique de la carte (par ex. #11/37.5508/-122.2895
  • Leaflet.label permet d’ajouter des étiquettes de texte aux marqueurs sur la carte et couches vectorielles
  • Leaflet.markercluster permet de manipuler des groupes de marqueurs de façon performante avec des animations Javascript

Le script charge ensuite les données de chaque couche, ces données ont été transformées au format Geojson par l’extension installée sous QGIS, ce sont ces données qui sont les sources de notre prototype par défaut.

Si l’on considère la couche qui comporte les informations départementales, le script va faire appel à la couche au format Géojson

<script src="data/exp_Departement.js"></script>

qui a pour contenu

var exp_Departement = {
   "type": "FeatureCollection",
   "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
   "features": [
 { "type": "Feature", "properties": { "color_qgis2leaf": 'none', "border_color_qgis2leaf": '#00ff00', "radius_qgis2leaf": 1.65, "transp_qgis2leaf": 0.75, "transp_fill_qgis2leaf": 1.0,   "id_geofla": 1, "code_dept": "01", "nom_dept": "AIN", "code_chf": "053", "nom_chf": "BOURG-EN-BRESSE", "x_chf_lieu": 8717, "y_chf_lieu": 65696, "x_centroid": 8814, "y_centroid": 65582, "code_reg": "82", "nom_region": "RHONE-ALPES" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 5.831226413621037, 45.938459578293219 ], [ 5.82212102507634, 45.930135953523816 ], [ 5.829154566881813, 45.917237123525148 ], ... } },
                   ...

GeoJson est un format qui permet d’encoder des structures de données géographiques, nous ne nous attarderons pas sur le format, notre objectif est de remplacer ces données statiques par des données fournies dynamiquement par notre serveur cartographique.

Une fois les données chargées, le programme passe au travail d’initialisation et de rendu de la carte.

Initialisation et rendu de la carte

var map = L.map('map', {zoomControl: true}).fitBounds([[36.1453640773, -8.27037589111], [54.0946095281, 6.49685600911]]);

Le script commence par créer un objet carte L.map, il passe en paramètre l’identifiant de l’élément div qui va contenir la représentation ainsi qu’un tableau d’options. L’option zoomControl permet d’indiquer que le contrôle de gestion de l’agrandissement sera intégré à la carte.

La méthode fitBounds() est chaînée à la méthode de création, elle définit la vue de carte qui contient les limites géographiques données avec le niveau de zoom maximum possible.
Les limites sont exprimées en couples de points, latitude et longitude.

var hash = new L.Hash(map);

Cette instruction permet d’utiliser l’extension leaflet-hash pour ajouter dynamiquement une ancre nommée à l’url fonction de la vue en cours qui permet de partager facilement l’état de la vue sur la carte. Elle peut être ignorée.

Le script va ensuite construire la liste des couches à afficher, deux instructions se trouvent dans le code, l’une utilise un featureGroup, l’autre un LayerGroup. Seule la première est exploitée.

var feature_group = new L.featureGroup([]);

Va déclarer un groupe qui permet de manipuler plusieurs couche à la fois au sein d’un seul élément. Ce groupe gère les événements souris.

La première couche sert à la représentation de la carte de base, le fond de carte.

            var basemap_0 = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
                attribution: additional_attrib + '&copy; <a href="http://openstreetmap.org">OpenStreetMap</a> contributors,<a href="http://creativecommons.org/licenses/by-sa/2.0/">CC-BY-SA</a>'});
            basemap_0.addTo(map);

L.tileLayer indique de charger une source de données depuis un serveur de tuiles (images pré-rendues).
La source de données est configurée pour interroger les serveurs de tuiles d’OpenStreetMap.
Le paramètre attribution est une chaîne de caractère représentatif qui sera affiché dans le contrôle qui gère les attributions
Ces attributions font parties des conditions d’utilisation d’OSM.

La couche est ajoutée directement à la carte.

Chaque couche de données est alors crée l’une à la suite de l’autre avec une source de données, un style et d’autres paramètres puis elle est ajoutée au groupe… qui au final n’est pas non plus utilisé ! L’extension sous QGIS manque d’aboutissement, et on doit faire abstraction de beaucoup de code généré. A chaque ajout de couche, la couche est temporisée dans un tableau, toutes les couches de la carte sont retirées et réinsérées dans l’ordre du tableau… passons, nous ne conserverons que les instructions d’ajout de couche à la carte.

Si l’on considère la couche des départements.

            var exp_DepartementJSON = new L.geoJson(exp_Departement, {
                onEachFeature: pop_Departement,
                style: function (feature) {
                    return {color: feature.properties.border_color_qgis2leaf,
                        fillColor: feature.properties.color_qgis2leaf,
                        weight: feature.properties.radius_qgis2leaf,
                        opacity: feature.properties.transp_qgis2leaf,
                        fillOpacity: feature.properties.transp_qgis2leaf};
                }
            });
...
map.addLayer(exp_DepartementJSON);

Ces instructions permettent de créer une couche de données vecteur, rendue côté client. La source de la couche est au format GeoJSON, elle accepte un objet GeoJSON, exp_Departement, qui a été construit lors de l’export du projet QGIS et chargé en amont dans les ressources Javascript.

L.geoJson() étend la classe featureGroup qui permet de gérer les événements souris.

Le style de rendu de la couche vectorielle est déclaré lors de la construction, il est fonction des éléments de la couche vectorielle, il est donc appliqué individuellement à chaque objet de la couche. Par défaut le code récupère le style depuis les propriétés exportées de l’extension Qgis2leaf, par exemple dans le Geojson vous pouvez trouver une propriété border_color_qgis2leaf positionnée à la couleur #00ff00.

Je peux modifier le style pour impacter tous les éléments avec pour l’exemple une couleur de remplissage rouge et une opacité qui permet de superposer cette couche au fond de carte.

                style: function (feature) {
                    return {color: feature.properties.border_color_qgis2leaf,
                        fillColor: '#FF0000',
                        weight: feature.properties.radius_qgis2leaf,
                        opacity: feature.properties.transp_qgis2leaf,
                        fillOpacity: 0.2};
                }

Remarquez alors le popup qui se déclenche au clic souris sur un département.

L’option onEachFeature permet d’indiquer une fonction qui sera appelée à la création d’une caractéristique vectorielle (chaque objet de la couche déclaré dans le format Geojson), ce qui permet d’attacher des événements aux caractéristiques vectorielles. Dans le code généré, une fonction est définie pour chaque couche, ici pop_Departement. Ce code ne fonctionne sur mon navigateur que si je modifie les styles pour l’affichage des départements. Il permet de manipuler chaque élément, par exemple d’afficher la propriété qui correspond au nom du département au clic sur le département sur la carte.

            function pop_Departement(feature, layer) {
                var popupContent = 'Departement ' + feature.properties['code_dept'];
                layer.bindPopup(popupContent);
            }

La fonction prend en argument les données géométriques et la couche vectorielle. Il est alors possible de récupérer les propriétés de l’objet, la méthode bindPopup() attache un événement au clic qui affichera le contenu html en paramètre dans un popup.

Ajout de contrôles à la carte

Un contrôle L.control( est un widget représenté à l’écran. Un contrôle peut être seulement informatif ou bien peut permettre d’interagir avec l’utilisateur. Nous avons vu lors de la création qu’un contrôle qui permettait de manipuler le niveau d’agrandissement était ajouté lors de la déclaration. Si vous jetez un œil à la carte vous verrez en fait quatre autres contrôles : le titre, l’échelle, les attributions et le contrôle qui permet de sélectionner les couches à afficher ou masquer.

Ces contrôles sont initialisés et ajoutés manuellement à la carte.

Contrôle pour affichage du titre

            var title = new L.Control();
            title.onAdd = function (map) {
                this._div = L.DomUtil.create('div', 'info');
                this.update();
                return this._div;
            };
            title.update = function () {
                this._div.innerHTML = '<h2>Etablissements scolaires de Gironde a 50 m parcelle cultivee</h2>Test client web leaflet'
            };
            title.addTo(map);

Le titre n’est pas un contrôle prédéfini, c’est un contrôle basique qui va implémenter l’interface IControl, qui lui donne accès à la méthode onAdd(). L’aspect visuel du contrôle est construit à l’ajout du contrôle à la carte, le but de la méthode d’ajout est de créer tous les éléments du DOM nécessaire pour le contrôle et d’ajouter des écouteurs d’événements, elle renvoie l’élément contenant le contrôle.

L’instruction L.DomUtil.create(‘div’, ‘info’) va créer dans le DOM un élément DIV de classe info.

L’élément est ensuite manipulé pour insérer un titre html.

Le système est souple et accorde beaucoup de liberté, la position du contrôle est quand à elle fournie en option à la création, à l’ajout à la carte, ou encore précisée post création.

La position est paramétrable dans les angles (haut droit, bas droit, haut gauche, bas gauche). Les contrôles sont rendus par des styles css qui permettent de gérer les marges par exemple, remarquez que le contrôle d’affichage du titre et le contrôle de visualisation des couches sont tous deux dans le coin haut droit, et empilés l’un au dessous de l’autre.

Contrôle pour affichage de l’échelle

L.control.scale({options: {position: 'bottomleft', maxWidth: 100, metric: true, imperial: false, updateWhenIdle: false}}).addTo(map);

L.control.scale() est un contrôle prédéfini qui permet l’affichage d’une échelle.

Contrôle pour visualisation des couches

L.control.layers(baseMaps, {"etabparcelinters50m": exp_etabparcelinters50mJSON, "etabssel": exp_etabsselJSON, "etabstamponselec": exp_etabstamponselecJSON, "DepartementWFS": exp_DepartementWFSJSON, "Departement": exp_DepartementJSON}, {collapsed: false}).addTo(map);

L.control.layers() est un contrôle prédéfini qui permet de sélectionner les couches représentées ou pas. Ce contrôle possède une option collapsed qui permet de représenter le contrôle sous forme d’icône cliquable, le contrôle s’ouvre ou se referme au clic.

Contrôle pour affichage des attributions

Le contrôle est créé et positionné par défaut sur la carte à moins de demander explicitement le contraire lors de la création de carte grâce au paramètre attributionControl positionné à faux. Les attributions sont récupérées des options de chaque couche. (La valeur pour la couche OSM est renseignée à la création)

L.control.attribution() permet de créer le contrôle manuellement, afin par exemple de modifier sa position.

Adaptation du code pour une source de données par serveur cartographique

Nous avons a disposition pour le projet deux services qui tournent sous un serveur cartographique QGIS Server.

Les capacités WFS et WMS des services peuvent être interrogées par les url suivantes

http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetCapabilities
http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

C’est ce dernier service que l’on va exploiter en premier.

Appel à un service WMS QGIS Server sous Leaflet

Le nom de la couche des départements nous intéresse particulièrement, on peut le récupérer depuis le fichier xml généré par la requête aux capacités du serveur.

<Layer queryable="1">
<Name>Departement</Name>
<Title>Departement</Title>
<CRS>EPSG:4326</CRS>
<EX_GeographicBoundingBox>
<westBoundLongitude>-5.21251</westBoundLongitude>
<eastBoundLongitude>9.63332</eastBoundLongitude>
<southBoundLatitude>41.3141</southBoundLatitude>
<northBoundLatitude>51.1376</northBoundLatitude>
</EX_GeographicBoundingBox>
<BoundingBox CRS="EPSG:4326" maxx="51.1376" minx="41.3141" maxy="9.63332" miny="-5.21251"/>
<Style>
<Name>default</Name>
<Title>default</Title>
<LegendURL>
<Format>image/png</Format>
<OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetLegendGraphic&LAYER=Departement&FORMAT=image/png&STYLE=default&SLD_VERSION=1.1.0"/>
</LegendURL>
</Style>
</Layer>

Leaflet est en mesure de charger des sources de données depuis un serveur cartographique via le service WMS, il utilise à cette fin TileLayer.WMS.

Nous pouvons remplacer la source de données Geojson par notre source WMS, la déclaration de la source se fait ainsi :

            var wmsDepts = L.tileLayer.wms("http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi", {
                layers: 'Departement',
                format: 'image/png',
                transparent: true,
            });

Le serveur cartographique retourne une image, de ce fait on perd les possibilités offertes par la couche vectorielle Geojson. L’avantage d’une image est souvent sa légèreté sur le réseau par rapport à un fichier de données, mais le format vectoriel nous permet d’interagir en javascript avec les objets vectoriels, on peut par exemple sélectionner un objet et modifier son style, supprimer ou ajouter un objet.

Interagir avec les objets peut permettre de créer des cartes plus interactives, en créant des composants de formulaire en widget par exemple on peut imaginer de représenter dynamiquement une carte choroplèthe ou des symboles proportionnels en fonction d’un paramètre sélectionné et d’une requête asynchrone qui va chercher les valeurs d’indicateurs à représenter.

Appel à un service WFS QGIS Server sous Leaflet

Si l’on regarde attentivement le code généré par l’extension sous QGIS, la solution est déjà implémentée mais ne fonctionne pas, le code tente non pas d’accéder au service WMS mais bien au service WFS paramétré sous QGIS. Il va effectuer une requête asynchrone en javascript sur le service WFS (de type GetFeature), le résultat est exploité comme couche vectorielle.

Comment s’y prend-il ?

            var DepartementWFSURL = 'http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=Departement&SRSNAME=EPSG:4326&outputFormat=text%2Fjavascript&format_options=callback%3AgetDepartementWFSJson';
            DepartementWFSURL = DepartementWFSURL.replace(/SRSNAME\=EPSG\:\d+/, 'SRSNAME=EPSG:4326');

Premièrement l’url d’appel au service WFS est construite. Cette url comporte le paramètre format_options qui ne fait pas parti du standard WFS, mais est interprété sur un serveur cartographique Geoserver. Ce paramètre permet au serveur d’encapsuler les données au format JSON dans une fonction javascript, les données retournées sont donc indiquées au format de sortie text-javascript. Ce mécanisme, nommé sous l’appellation JSONP, permet de contourner les appels de requêtes sur domaines croisés.

Si l’on effectue une requête sur notre serveur afin de récupérer des données JSON depuis une page html côté client (en Ajax donc), le navigateur va interdire l’opération si le serveur n’autorise pas explicitement les accès en origine croisée par l’ajout d’une entête de contrôle d’accès.

Notre problème est double. Premièrement QGIS Server possède les capacités pour servir du GeoJSON, qui est simplement du JSON avec une convention, par contre il ne propose pas le format de sortie JSONP qui encapsulerait simplement le GeoJSON dans une fonction Javascript. Il nous resterait comme option de configurer le serveur pour accepter les requêtes d’origine croisées (CORS), cependant le serveur QGIS fonctionne en mode CGI et Apache ne peut pas à ma connaissance ajouter une entête d’autorisation sur la sortie produite par QGIS Server.

La solution lorsque l’on ne peut ni modifier le format de sortie ni configurer le serveur pour ajouter des autorisations est d’utiliser un proxy. Très simplement un script sera hébergé sur le même serveur, il acceptera en paramètres un nom de fonction qui servira à l’encapsulation javascript du format JSON, et l’url qu’il sera en mesure d’interroger depuis le même domaine.

Corrigée pour le serveur QGIS, l’url qui permet de récupérer les objets de la couche des départements est la suivante :

var DepartementWFSURL = 'http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WFS&VERSION=1.0.0&REQUEST=GetFeature&TYPENAME=Departement&SRSNAME=EPSG:4326&outputFormat=GeoJSON';

J’ai créé un simple script en langage PHP qui servira de proxy. Vous pourrez probablement trouver des solutions plus complètes…

<?php
if(!isset($_GET['callback']) || !isset($_GET['url'])) {
  header('status: 400 Bad Request', true, 400);
  exit;
}
extract(parse_url($_GET['url']));       

if (! ( in_array($host, array($_SERVER['HTTP_HOST'], $_SERVER['SERVER_ADDR']) )
      || ( basename($path) != 'qgis_mapserv.fcgi' ) )) {
  header('status: 400 Bad Request', true, 400);
  exit;
}
if (!$content = file_get_contents(filter_var($_GET['url'], FILTER_SANITIZE_URL))) {
  header('status: 400 Bad Request', true, 400);
  exit;
}
header('content-type: application/javascript; charset=utf-8');
header("access-control-allow-origin: *");
echo filter_var($_GET['callback'], FILTER_SANITIZE_ENCODED), '(', $content, ');';

Ce script nommé proxy.php est placé à la racine du répertoire web. Il doit être appelé avec deux paramètres, callback et url.

Le paramètre url doit être encodé…

http://io.gchatelier.fr/proxy.php?callback=getDepartementJson&url=...

Le script effectue quelques contrôles, premièrement il n’aboutit pas si les paramètres ne sont pas renseignés. Il vérifie ensuite que l’url à interroger est sur le même domaine que le serveur, puis vérifie que l’on cherche bien à interroger un script CGI qui correspond à QGIS Server. Le script filtre ensuite les paramètres et exécute la requête. Le contenu de la requête est retourné encapsulé dans la fonction de rappel Javascript avec les bons entêtes : premièrement l’entête d’application javascript (et non pas l’entête d’application JSON), puis l’entête de contrôle d’accès qui autorise le proxy à être interrogé depuis un autre domaine.

Bien… ce problème réglé le code peut fonctionner comme attendu.

            var proxy = 'http://io.gchatelier.fr/proxy.php?callback=getDepartementJson&url=';
            var proxyURL = proxy + encodeURIComponent(DepartementWFSURL);

Construit l’url d’appel à notre proxy en précisant la fonction de rappel et l’url d’origine du service.

            var exp_DepartementWFSJSON = L.geoJson(null, {
                style: function (feature) {
                    return {color: '#00ff00',
                        fillColor: '#ff0000',
                        weight: 1.65,
                        opacity: 0.75,
                        fillOpacity: 0.10};
                },
                onEachFeature: function (feature, layer) {
                    var popupContent = 'Departement ' + feature.properties['code_dept'];
                    layer.bindPopup(popupContent);
                }
            });

Le code commence par créer une couche GeoJSON vide. Le mécanisme est le même que nous avons rencontré précédemment, avec application d’un style (sur mon navigateur le style de remplissage doit être renseigné pour que le popup fonctionne), et application d’une fonction javascript à exécuter sur chaque objet de la couche vectorielle. Comme précédemment on créé un popup au clic souris sur le département pour afficher son code.

Cette couche vide sera remplie ultérieurement par un appel javascript asynchrone qui va récupérer les objets de la couche vectorielle depuis le service WFS.

            var DepartementWFSajax = $.ajax({
                url: proxyURL,
                dataType: 'jsonp',
                jsonpCallback: 'getDepartementJson',
                success: function (response) {
                    L.geoJson(response, {
                        onEachFeature: function (feature, layer) {
                            exp_DepartementWFSJSON.addData(feature)
                        }
                    });
                }
            });

            exp_DepartementWFSJSON.addTo(map);

L’appel Ajax accepte ces paramètres :

  • url : url vers notre proxy qui interrogera le serveur WFS
  • dataType : notre proxy retourne du javascript au format JSONP, on doit préciser le type de données car les entêtes retournées sont application javascript
  • jsonpCallback : le nom de la fonction de rappel dans le fichier JSONP retournée par le proxy
  • success : En cas de réussite une fonction anonyme est exécutée avec la réponse en paramètre, notre GeoJSON… la réponse est parcourue dans une couche L.geoJson qui duplique chaque objet dans la couche vide

Le développeur de l’extension d’export sous QGIS a fait le choix de dupliquer la couche en cas de réussite ce qui permet d’initialiser une couche vectorielle y compris au cas où la requête Ajax échoue.

Nous avons fait le tour des codes sources générés sous une extension QGIS pour la librairie Leaflet. Vous pouvez également vous référer au même document adapté pour la librairie OpenLayers3.

Nous sommes désormais en mesure de construire des applications interactives pour le web depuis des services cartographiques, ce sera l’objet du prochain document de cette série qui détaillera la création d’une carte interactive sous la librairie Leaflet.

Utilisation des services WMS et WFS sous un client web Leaflet ou OpenLayers 3

Utilisation des services sous un client web

Nous avons vu comment créer une carte sous le logiciel client QGIS et comment publier le projet sous le serveur cartographique QGIS server au travers de deux types de services, WMS et WFS. Nous sommes donc en mesure de servir des données cartographiques au travers du Web. Ce n’est que la face cachée de l’iceberg, nous allons exploiter ces services pour construire une application de rendu de carte dynamique dans un navigateur web.

Les solutions clientes qui permettent le rendu de cartes interactives sur un navigateur sont des librairies en langage Javascript. Nous couvrirons deux de ces librairies : premièrement la librairie Leaflet qui a été conçue après l’apparition du HTML5 dans un souci de simplicité et d’efficacité, puis on s’attardera sur la référence dans le domaine, la librairie OpenLayers. Pour cette dernière j’emploierai la version 3 qui est une refonte moderne de la version 2, désormais obsolète pour les futur développements.

QGIS possède des extensions pour ces deux librairies qui permettent de générer un squelette de projet web très simplement à partir du projet en cours. Nous allons utiliser ces extensions qui nous facilitent le travail pour initialiser avec les bons paramètres les couches du projet à charger. Ces extensions permettent de générer dans un fichier au format Json les données de couches exploitées par la librairie javascript, alors que l’on souhaite interroger le serveur cartographique. Nous testerons la capacité de ces extensions à générer du code d’appel aux services WFS puis nous verrons en détail le code source nécessaire pour interroger un serveur distant.

Rappels sur le projet de test

Le projet de test est celui que nous avons publié sur le serveur cartographique. Il comporte cinq couches destinées à mettre en évidence les intersections des collèges et lycées avec des parcelles cultivées dans un rayon de 50 mètres, sur le département de la Gironde.

Les couches Departement et cultures_33 servent à représenter les limites administratives, et limites de parcelles. Elles sont récupérées depuis la base de données. Les autres couches sont des couches projet enregistrées au format fichier Shapefile, elles servent à matérialiser les établissements que l’on a identifiés : la couche etabs_sel permet d’afficher un point aux coordonnées de latitude et longitude de l’établissement, la couche etabs_tampon_selec sert à délimiter le périmètre autour de l’établissement et finalement la couche etab_parcel_inters_50m met en évidence l’intersection entre la parcelle de culture et la zone tampon de 50 mètres autour de l’établissement.

Préparation du projet de test

Dupliquons dans le projet la couche des départements chargée dans la base locale avec la même couche mais chargée cette fois depuis le serveur cartographique.
Cette manipulation permettra de tester la capacité des extensions QGIS à générer le code client pour l’exploitation du service.

Ouvrez le projet QGIS que vous avez publié sur le serveur cartographique et enregistrez le sous un autre nom
(ex. C:\SIG\projet_1\projet_1.qgs vers C:\SIG\projet_1\projet_1_distant.qgs)

Ajoutez la couche WFS des départements (Depuis l’icône à gauche ou bien Couche > Ajouter une couche WFS)

Vous devriez avoir désormais deux couches de départements… il vous suffit de copier le style de la première et de l’appliquer à la seconde.

Faites un clic droit dans l’explorateur de couche sur la couche initiale, sélectionnez l’entrée de menu Copier le style.
Procédez de même sur la couche de destination pour Coller le style.

Modifiez le nom de la couche de destination puis enregistrez !

Préparation du projet sous QGIS pour le client Leaflet

Sous QGIS, ouvrez le gestionnaire d’extension depuis le menu principal
Extension > Installer / Gérer les extensions
Recherchez et installez l’extension qgis2leaf
Ouvrez le projet QGIS que vous avez publié sur le serveur cartographique.
(ex. C:\SIG\projet_1\projet_1.qgs)

Depuis le menu principal vous pouvez alors créer une carte web
Internet > qgis2leaf > Export a QGIS project to a working leaflet webmap

Je n’ai renseigné qu’un minimum d’options afin dans un premier temps de générer le moins de code possible :

  • Appuyer sur le bouton Get Layers permet de charger vos couches dans l’extension. Vous pouvez sélectionner les couches à exporter, ici j’ai sélectionné toutes les couches
  • Le champ Frame width / height permet de déterminer les dimensions de la carte sur la page web, j’ai sélectionné le champ Full screen pour utiliser tout l’écran
  • Le champ Extent permet de définir l’étendue géographique initiale à l’arrivée sur l’application, par défaut elle est positionnée à l’étendue du canvas en cours, ici je me suis positionné sur une vue globale centrée sur la Gironde
  • Le champ Basemaps permet de faire une sélection multiple de cartes à importer pour le fond de carte, pour l’habillage j’utilise le fond de carte OpenStreetMap OSM Standard

Le projet est enregistré et ouvert dans un navigateur web, le résultat est fonctionnel on peut masquer des couches et zoomer. Par contre la couche des cultures est affichée quelle que soit l’étendue géographique, alors que l’on avait paramétré dans le projet sa visibilité. Le chargement très long de la page est lié à la quantité de données de cette couche chargée intégralement et sera réduit lorsqu’on chargera la couche depuis le service web en fonction de l’étendue géographique. Notez également que la couche WFS des départements n’est pas rendue…

Pour terminer avec les problèmes relevés : l’étendue initiale du canevas n’est pas respectée.

Ce premier rendu est prometteur, nous rentrerons dans un second temps dans le code qui a été généré, pour le moment voyons en parallèle comment arriver aux mêmes résultats avec une autre librairie.

Préparation du projet sous QGIS pour le client OpenLayer 3

Sous QGIS, ouvrez le gestionnaire d’extension depuis le menu principal
Extension > Installer / Gérer les extensions
Dans l’onglet Paramètres cochez la case pour afficher les extensions expérimentales.
Recherchez et installez l’extension Export to OpenLayers 3
Ouvrez le projet QGIS préparé précédemment. Vous pouvez alors créer une carte web OpenLayers3
Internet > Export to Open Layers > Create OpenLayers map

J’ai laissé les paramètres par défaut devant le peu de choix et l’absence d’aide à l’utilisation. L’extension s’exécute, cela requiert du temps à nouveau en raison de la couche des cultures qui représente 20 Mo de données au format Json. La carte possède les mêmes qualités et défauts que la carte générée pour la librairie Leaflet, c’est à dire que les cultures sont représentées indépendamment de l’échelle, le service WFS n’est pas non plus exploité car les données ont été converties au format Json.

Nous avons vu comment générer les squelettes des codes sources qui permettent d’utiliser des clients web pour l’affichage cartographique. Pour chaque solution nous allons explorer plus en détail le code afin d’expliquer le fonctionnement des librairies et comment nous pouvons remplacer les données générées en dur par des données issues des services que l’on a mis en place.

Installation de QGIS Server et publication d’un projet QGIS

QGIS Server

Introduction

QGIS Server est une application qui va agir en service et répondre à des requêtes http qui suivent un standard conçu pour les échanges de données géospatiales sur le web.

Il s’installe en application FastCGI/CGI sur un serveur web et va générer des réponses aux formats standards qui permettrons d’échanger des données, de visualiser et publier des cartes sur Internet à partir de données géographiques.

Ce serveur est conçu à partir des librairies du projet QGIS, il possède le même moteur de rendu cartographique que l’application QGIS Desktop et permet de publier très simplement un projet QGIS Desktop vers le Web.

Le serveur cartographique possède 2 services nommés WMS et WFS que l’on utilisera au travers de clients web (Leaflet et Openlayers 3)

Serveur WMS

QGIS Server implémente le standard WMS (Web Map Service) publié par l’Open Geospatial Consortium (OGC). WMS est un service web mis en oeuvre au travers d’une architecture REST, il est possible de l’interroger par des requêtes HTTP. Le serveur va répondre aux requêtes, l’objet du service WMS est simplement de renvoyer une image : le serveur analyse la requête, détermine les paramètres, et en fonction, il va construire une carte raster à partir des sources de données dont il dispose.

Les opérations basiques du WMS sont :

  • GetCapabilities : Retourne les métadonnées sur les capacités du service. (Les paramètres de requête acceptés par le serveur)
  • GetMap : Retourne une carte (au format d’image) avec des paramètres définis, dimensions et informations
  • GetFeatureInfo : Retourne des informations sur les entités d’une carte

Serveur WFS et WFS-T

QGIS Server implémente également le standard WFS (Web Feature Service) publié par l’OGC. Contrairement au WMS, le WFS n’implique pas le rendu d’une carte par le serveur, il s’agit en fait d’une requête de récupération d’informations sous forme vectorielle. Les entités d’une couche sous QGIS peuvent être communiquées sous forme xml au client qui en fait la demande. Le client du service peut alors utiliser le service WFS pour effectuer par exemple un rendu de carte.

WFS-T est une version transactionnelle qui autorise la création, suppression, et mise à jour d’entités.

Référencez vous aux documents de standards pour les paramètres de construction de requêtes http

Et quoi d’autre ?

Eh bien à vous de découvrir… nous avons presque tout ce qu’il faut pour construire une application cartographique, les autres capacités du serveur ont peu d’importance dans ce contexte.

Nous verrons à partir d’un client web comment interroger les serveurs de tuile d’OpenStreetMap pour le rendu d’un fond de carte, QGIS Server n’a pas les capacités de serveur de tuiles à l’heure actuelle. En fin de projet je testerai l’installation des données OSM sur un serveur qui ne sera dédié qu’à cette opération. Il pourra être utilisé en remplacement des serveur OSM afin de ne pas pénaliser ceux ci.

Installation

Sous Windows

Vous pouvez installer QGIS Server en local sur votre poste de travail Windows, vous pouvez utiliser l’outil d’installation OSGeo4W fourni par la Fondation Geospatiale Open Source (OSGeo)

L’installateur propose une liste de paquets logiciels qui permet d’installer les logiciels et leurs dépendances. Le serveur QGIS se situe sous l’arborescence Web.

L’installation embarque une version du serveur web Apache 2 qui est installée comme un service, lancé au démarrage. Vous pouvez gérer le serveur Apache depuis une entrée du menu Démarrez.

Après installation vous devrez renseigner le port d’écoute d’Apache, particulièrement si vous possédez déjà des serveurs par défaut sur le port 80.

Il suffit d’éditer le fichier C:\OSGeo4W\apache\conf et remplacer la ligne

Listen @apache_port_number@

par votre port d’écoute. Par ex. sur le port 80

Listen 80

Redémarrez Apache depuis l’entrée de menu
Tous les programmes > OSGeo4W > Apache

Sous Debian Wheezy

Debian Wheezy possède des paquets compilés disponible via l’outil d’installation, ce qui rend la tâche en théorie triviale, cependant nous verrons plus loin que le paquet par défaut propose une ancienne version du serveur qui ne possède pas le support WFS.

La commande installera également Apache 2 avec un module capable de gérer le FastCGI.

aptitude install qgis-mapserver libapache2-mod-fcgid

Autoriser le module CGI

a2enmod cgid

Redémarrer Apache

service apache2 restart

Et c’est tout ? Oui. Sous Debian, et pour la version du serveur qui ne propose que le service WMS.

Continuez la lecture ! Vous trouverez plus loin les instructions pour installer la dernière version stable afin de profiter du support WFS

J’ai moi même procédé en 2 étapes et je préfère conserver une trace du cheminement complet afin d’aider les personnes qui pourraient se poser les mêmes questions.

Le script qgis_mapserv.fcgi est placé sous le répertoire /usr/lib/cgi-bin/, c’est lui que l’on va interroger.

Tester l’installation

Le serveur tourne, il doit répondre aux requêtes selon des standards précis. On peut interroger le serveur et lancer une requête http GetCapabilities sur les services WMS et WFS pour qu’il nous informe de ce qu’il est capable de faire.

Par défaut sous Windows

http://localhost/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities
http://localhost/qgis/qgis_mapserv.fcgi.exe?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetCapabilities

Par défaut sous Debian

http://92.222.27.205/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities
http://92.222.27.205/cgi-bin/qgis_mapserv.fcgi?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetCapabilities

Nous allons voir comment ajouter un premier projet QGIS sur le serveur.

Mettre en ligne un projet QGIS Server

Préparation du projet sous QGIS

L’objet de ce premier projet sera de vérifier le fonctionnement du serveur, on se contentera d’afficher sur une carte les départements et les établissements scolaires à moins de 50 mètres d’une parcelle agricole sur le département 33

Ouvrez un nouveau projet sous QGIS. (Par défaut votre SCR doit être EPSG:4326)
Chargez les couches PostGIS des établissements, des départements, et de la parcelles pour le département 33 (GIRONDE) :

Couche > Ajouter une couche PostGIS

Le style de rendu des départements doit surtout mettre en avant le contour de polygone pour l’information de limite administrative.

Le contenu du polygone n’apporte pas d’information dans ce cas aussi nous allons le rendre opaque et accentuer les contours.

Ouvrez les propriétés de la couche Departement. (Il suffit de faire un double clic sur le nom de la couche dans l’explorateur de couche)
Dans la fenêtre de propriétés, ouvrez l’onglet Style.

Paramétrez le style de remplissage sur Bordure : Ligne simple pour le type de symbole.

Vous pouvez également modifier la couleur de remplissage et la transparence de la couche pour améliorer le rendu.

A ce stade vous pouvez enregistrer le projet (ex. C:\SIG\projet_1\projet_1.qgs)

La distance de 50 mètres autour d’un établissement peut se matérialiser approximativement par une zone tampon autour du point. Une zone tampon constitue une délimitation autour d’un objet géographique qui épouse la forme de l’objet. Autour d’un point une zone tampon aboutira plus ou moins à un cercle ayant pour centre l’objet point et pour rayon la distance tampon. La création d’une zone tampon permettra de déterminer s’il y a une intersection du tampon avec les objets de type polygone qui représentent les parcelles.

Ouvrez la fenêtre de création de tampons :
Vecteur > Outils de géotraitement > Tampon(s)

Renseignez le formulaire :

  • Couche vectorielle de saisie : il s’agit de la couche autour de laquelle seront crées les zones tampon
  • Segments pour l’approximation : plus la valeur est élevée plus les contours de la zone tampon seront arrondis
  • Distance tampon : Le système de coordonnées de référence du projet est WGS 84, l’unité du canevas est de ce fait exprimée par défaut en degrés décimaux ce qui n’est pas pratique pour matérialiser une distance en mètres ! Une valeur de 0.0005 donnera un résultat acceptable pour ce simple prototype. C’est très dommage que QGIS n’indique pas l’unité utilisée à cet endroit
  • Fichier de sortie : Choisissez l’emplacement de la couche qui sera créée, par exemple sous le répertoire projet C:/SIG/projet_1/etab_tampon_50m.shp
  • Ajouter le résultat au canevas

Si l’on utilise l’outil de mesure inclus dans QGIS on peut se rendre compte que le rayon de la zone tampon varie autour de 50 mètres (entre 40 et 56 mètres), c’est plus flagrant si l’on diminue le nombre de segments pour l’approximation dans la fenêtre de paramétrage de création du tampon.

L’intersection est l’opération qui va permettre de créer une couche avec les parties communes à la couche tampon autour des établissements et la couche de parcelles.

Ouvrez la fenêtre de création d’intersection :
Vecteur > Outils de géotraitement > Intersection
Dans le formulaire qui s’ouvre alors, sélectionnez la couche des tampon et celle des parcelles de la Gironde puis ajoutez le résultat au canevas

N’essayez pas sur les parcelles pour la France entière, vous perdriez certainement votre temps…

Il reste alors à distinguer les établissements pour lesquels il existe une intersection. On peut pour cela ajouter une jointure sur le numéro d’établissement entre la couche des établissements et la couche des intersections puis à partir des attributs appliquer un filtre de sélection qui va restreindre les éléments sélectionnés sur la base de la base de l’existence des attributs de la table jointe.

Créer une jointure dans les propriétés de la couche établissement


Procéder à la sélection dans la table des attributs à l’aide d’un filtre avancé sous forme d’expression

NOT  "etab_parcel_inters_50m_appellatio"  IS NULL

Il suffit alors d’enregistrer la sélection en sauvegardant la couche établissement sous un autre nom et en précisant de n’enregistrer que les entités sélectionnées.

A ce stade j’ai supprimé la couche d’établissements du projet et la couche de tampons, puis appliqué un style à la nouvelle couche d’établissements à moins de 50 mètres d’une parcelle en Gironde.

Dernier point, j’ai rendu la visibilité des parcelles dépendante de l’échelle afin de n’afficher les parcelles que lorsque l’utilisateur aura suffisamment zoomé sur une école. Cela peut se faire à partir des propriétés générales de la couche.

Ce projet est presque prêt à être publié, il reste à paramétrer les propriétés du serveur dans les propriétés du projet.

Publication du projet sous QGIS Server

Paramétrer le projet pour publication

Ouvrez les propriétés du projet (CTRL + MAJ + P)
Projet > Propriétés du projet

Sous l’onglet Serveur OWS vous trouverez quatre rubriques, les plus importantes pour nous étant les capacités WMS et WFS

Sous la rubrique Informations générales du service vous pouvez instruire des informations qui seront communiquées par le service, c’est utile si vous destinez les services à des clients web externes à votre organisation. Vous pouvez renseigner ici à minima les informations pour vous contacter.

Faites un clic droit sur la couche Departement, et sélectionnez l’entrée de menu Zoomer sur la couche

L’emprise du projet est adaptée en conséquence, nous avons alors une vue globale de la France métropolitaine, ce sera notre vue par défaut.

Sous la rubrique Capacités WMS, cochez la case Emprise annoncée, et cliquez sur le bouton Utiliser l’emprise actuelle du canevas

Restreignez les SCR au SCR utilisé pour le projet (EPSG:4326), sans quoi le service listera tous les SCR lors d’une requête d’information sur ses capacités.

Sous la rubrique Capacités WFS vous pouvez sélectionner les couches que vous souhaitez publier en WFS, les colonnes Mise à jour, Insérer et Effacer permettent de sélectionner finement les opérations transactionnelles à autoriser (WFS-T).

Cliquez sur le bouton Sélectionner tout pour publier toutes les couches. Nous ne souhaitons pas publier en mode transactionnel, vous pouvez appliquer et enregistrer le projet

A ce stade le projet peut être mis en publication sous le serveur QGIS mais il comporte encore beaucoup d’informations qu’il n’est pas utile de transmettre sur une connexion web.

Pour la plupart des couches il est inutile de transmettre tous les attributs (c’est d’autant plus vrai que l’on a à chaque fois conservés tous les attributs lors des jointures), je supprime par défaut tous les attributs et n’autorise qu’un nombre restreint d’attributs qui pourraient utiles à l’application.

Faites un clic droit sur la couche Departement, et sélectionnez l’entrée de menu Propriétés

Sous l’onglet Champs QGIS dresse la liste des attributs d’une couche, avec pour dernières colonnes des cases à cocher qui permettent de publier ou pas l’attribut en service WMS et WFS.

Décochez les champs superflus, et procédez de même pour les autres couches.

Enregistrez le projet !

Cette fois ci, le projet est prêt ! Passons à la mise en route.

Prise en charge du projet par le serveur

Sous Windows

Pour ajouter le projet au serveur il suffit de créer un répertoire projet sous C:\OSGeo4W\apps\qgis\bin, copier l’exécutable CGI sous ce répertoire ainsi que le fichier projet *.qgs et les ressources qu’il utilise.

Lors de la création du projet j’ai enregistré le projet au même endroit que les couches créées, QGIS s’appuies sur un chemin relatif pour accéder aux ressources, il utilisera le chemin ./ pour accéder à chaque ressource, il vous incombe de respecter le chemin relatif vers vos données. Continuez la lecture, je vous montrerai plus loin comment modifier la source d’une couche pour un projet QGIS !

Si vous possédez une installation cygwin (Sinon c’est aussi rapide à la souris mais avec moins de valeur documentaire !):

cd C:/OSGeo4W/apps/qgis/bin
mkdir projet1
cp qgis_mapserv.fcgi.exe projet1/qgis_mapserv.fcgi.exe
cp C:/SIG/projet_1/* projet1

Vérifier le fonctionnement

WMS

http://localhost/qgis/projet1/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

Et pour ne pas vous décevoir, une carte… cette url est construite d’après les paramètres retournés par la requête GetCapabilities, à vous de fouiller pour trouver où.

http://localhost/qgis/projet1/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&CRS=EPSG:4326&REQUEST=GetMap&BBOX=-5.82882,41.1196,10.2496,51.3322&WIDTH=800&HEIGHT=800&LAYERS=Departement,etab_parcel_inters_50m,etabs_tampon_selec,etabs_sel&FORMAT=image/png&STYLE=

WFS

http://localhost/qgis/projet1/qgis_mapserv.fcgi.exe?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetCapabilities

Si cela vous amuse vous pouvez ouvrir un nouveau projet sous QGIS et tenter d’accèder au service :
Couche > Ajouter une couche WFS

Vous pouvez alors vous connecter, sélectionner toutes les couches à importer et vérifier le résultat.

Sous Debian

Pour les besoins du projet j’ai travaillé sur une connexion PostGIS locale, cela implique que le projet pourrait ne pas fonctionner si je le transfère sur un serveur.

Vous pouvez modifier les sources de données pour chaque couche en éditant le fichier projet *.qgs avec n’importe quel éditeur de texte, il s’agit simplement d’un format xml.

Il suffit de rechercher et adapter les lignes de ce type

<datasource>dbname='Pesticides' host=localhost port=5432 user='postgres' password='postgres' sslmode=disable key='id_geofla' srid=4326 type=MULTIPOLYGON table="public"."Departement" (geom) sql=</datasource>

ou encore

<datasource>./etab_parcel_inters_50m.shp</datasource>

Dans l’immédiat il me manque la table des parcelles sur la Gironde sur le PostgreSQL distant, l’utilitaire de sauvegarde et restauration de PostgreSQL pg_dump peut me permettre d’exporter cette table vers un fichier sql.

cd C:\SIG\sql
pg_dump --host=localhost --port=5432 --username=postgres --dbname=Pesticides --no-password --table=cultures_33 > db.sql

Sous le serveur distant, la restauration se fait par :

su - postgres
psql --dbname=Pesticides --file=db.sql

Une fois l’environnement dupliqué l’installation du projet est très simple.

Il suffit comme sous Windows de créer un répertoire et de copier sous ce répertoire le fichier projet *.qgs et les ressources. Un lien symbolique doit suffire pour donner accès au script CGI depuis le répertoire projet.

mv projet1 /usr/lib/cgi-bin/
cd /usr/lib/cgi-bin/projet1
ln -s ../qgis_mapserv.fcgi .
ln -s ../wms_metadata.xml .

Vérifier le fonctionnement

WMS

http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities
http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&CRS=EPSG:4326&REQUEST=GetMap&BBOX=41.1196,-5.82882,51.3322,10.2496&WIDTH=800&HEIGHT=800&LAYERS=Departement,etab_parcel_inters_50m,etabs_tampon_selec,etabs_sel&FORMAT=image/png
Pour construire cette dernière url j’ai du inverser l’axe x,y du paramètre BBOX, autrement le serveur retournait une image vide. Je ne trouves pas les informations pour paramétrer le serveur sous Debian, mais à priori cela correspond à la version 1.1.0 du service, le numéro de version semble tout simplement ignoré.
WFS : Attention problème ! Si vous testez le lien suivant le serveur répond avec les capacités WMS, le paramètre est simplement ignoré, la version du serveur est trop ancienne !
http://92.222.27.205/cgi-bin/projet1/qgis_mapserv.fcgi?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetCapabilities

Mais pourquoi ça ne fonctionne jamais du premier coup ?

Installer la dernière version stable de QGIS Server sur Debian Wheezy

En prérequis vous devrez avoir installé et configuré le paquet libapache2-mod-fcgid comme expliqué en amont.

Par défaut sous Debian Wheezy l’installeur ne propose pas la dernière version du serveur QGIS… en fait le paquet qgis-mapserver a été remplacé par le paquet qgis-server.

Retirez si besoin le paquet obsolète
apt-get remove qgis-mapserver

Il faut ensuite indiquer à Debian où récupérer les dernières versions stables du nouveau paquet.

Éditez le fichier de sources du gestionnaire de paquets
vi /etc/apt/sources.list.d/pgdg.list
Insérez les lignes suivantes : (Consultez la documentation QGIS pour votre version Debian)
deb     http://qgis.org/debian wheezy main
deb-src http://qgis.org/debian wheezy main
Récupérez la clef de sécurité et mettez à jour les informations dans le gestionnaire de paquets
gpg --keyserver keys.gnupg.net --recv-key DD45F6C3
gpg --export --armor DD45F6C3 | apt-key add -
apt-get update
Procédez à l’installation
apt-get install qgis-server
apt-get install build-essential python-dev

Cette dernière commande permet de régler des problèmes de fonctionnement, vous en aurez probablement également besoin.

Nous avions déjà tout paramétré, il ne vous reste qu’à reprendre à partir de l’étape de vérification du fonctionnement et vous pouvez désormais vérifier le chargement des couches WMS et WFS depuis QGIS.

En cas de problème

Les fichiers de log du script serveur qgis-server peuvent servir à résoudre des problèmes… encore faut-il savoir les trouver !

Sous Apache et avec le mode fcgid, vous pouvez éditer le fichier de configuration Apache pour insérer deux directives.

FcgidInitialEnv QGIS_SERVER_LOG_FILE /tmp/qgisserver.log
FcgidInitialEnv QGIS_SERVER_LOG_LEVEL 0

Remarquez qu’en dehors du répertoire temporaire je n’ai pas réussi à faire fonctionner ces directives.

Côté performances vous pouvez également paramétrer le mode fcgid.

vi /etc/apache2/mods-enabled/fcgid.conf

Ma configuration est la suivante :

<IfModule mod_fcgid.c>
  AddHandler    fcgid-script .fcgi
  FcgidConnectTimeout 5
  FcgidIOTimeout 240
  FcgidMaxProcesses 1500
  FcgidMaxProcessesPerClass 150
  FcgidMinProcessesPerClass 10
  FcgidMaxRequestInMem 655360
  FcgidMaxRequestLen 1310720
  DefaultInitEnv LD_LIBRARY_PATH /usr/lib:/usr/lib
</IfModule>

Avant PostgreSQL 9.3 il était aussi nécessaire de bien configurer le système pour des performances sur la base de données.

Pour référence les paramètres du kernel sont par défaut sous dimensionnés, ils peuvent être déterminés simplement par un script repris depuis ce Blog.

vi shmsetup
chmod o+x shmsetup
#!/bin/bash
page_size=`getconf PAGE_SIZE`
phys_pages=`getconf _PHYS_PAGES`
shmall=`expr $phys_pages / 2`
shmmax=`expr $shmall \* $page_size`
echo kernel.shmmax = $shmmax
echo kernel.shmall = $shmall
./shmsetup >> /etc/sysctl.conf
sysctl -p

Epilogue

Vous devez probablement vous dire que tout ce que nous retirons du serveur est pour l’instant bien austère et manque cruellement de dynamisme par rapport à ce que l’on peut voir ailleurs et vous avez raison mais l’on a rien sans rien, et bientôt vous serez comblé. Dans le prochain article nous verrons comment faire appel à ces services sous un client Web

Charger des données spatiales : alimenter PostgreSQL PostGIS

Nous avons récupéré les données spatiales qui vont servir à développer l’application cartographique puis installé la base de données et son extension spatiale qui vont servir à enregistrer les données.

Quand se pose la question de l’alimentation de la base, une multitude de solutions se profilent. On peut par exemple utiliser un outil logiciel en ligne de commande, une extension à un logiciel existant, QGIS ou Pg admin proposent de telles extensions, ou encore utiliser un ETL (Extract Transform Load) comme Geokettle afin d’automatiser le processus.

Nos jeux de données possèdent quelques particularités comme nous l’avons vu et je vais couvrir deux cas de figure pour répondre au besoin.

Gérer les tables au travers de QGIS Desktop

QGIS est une solution complète qui permet d’effectuer beaucoup d’opérations simplement à la souris et avec beaucoup de contrôle par rapport à un outil en ligne de commande. Je vais l’utiliser en priorité avec l’extension DB Manager.

Sur les anciennes version de QGIS l’extension SPIT permettait déjà l’import de fichiers au format Shapefile vers une base de données. Cette extension sera remplacée à terme par un système unifié, DB Manager que nous utiliserons pour nos besoins. (J’ai rencontré par ailleurs des problèmes d’import avec SPIT sur les dernières versions)

Import des données et coordonnées géographiques des établissements scolaires français

Ouvrez QGIS et procédez à l’ouverture du fichier de données (référez vous au document de préparation des données).

Le logiciel écarte automatiquement 120 enregistrements pour absence de géométrie, nous n’avons pas à faire l’opération manuellement.

Vous devriez alors vous retrouver avec une représentation des établissements par des points pour les établissements de métropole et d’outre mer.

Comme nous l’avons vu précédemment, le fichier possède plusieurs système de coordonnées pour les valeurs X et Y, et seules les valeurs pour la France métropolitaine et la Corse sont exprimées dans le système de coordonnées RGF93 / Lambert-93, les points pour les établissements en dehors de la métropole sont donc mal positionnés sur la carte.

Nous avons également vu qu’au niveau des attributs seul le code postal est un facteur discriminant, nous ne possédons pas d’information sur le département ou la commune.

Voyons comment éliminer simplement les établissements d’outre mer pour lesquels nous ne possédons pas de données parcellaires sur les cultures agricoles.

La liste des codes département à l’outre mer est :

971 Guadeloupe
972 Martinique
973 Guyane
974 Réunion
975 St Pierre et Miquelon
976 Mayotte
986 Wallis et Futuna
987 Polynésie Française
988 Nouvelle Calédonie

Le code postal pour la France d’outre mer commence par les 3 chiffres du code départements, on peut s’arrêter aux 2 premiers chiffres pour éliminer l’outre mer, la recherche d’un numéro 97 ou 98 en première position du code postal permet de sélectionner un enregistrement à éliminer.

Ouvrez la table d’attributs (clic droit sur la couche dans l’explorateur de couches).

Par défaut le tableau fait afficher tous les enregistrements. QGIS permet de faire une sélection d’enregistrements par filtre et d’enregistrer la sélection comme une nouvelle couche. C’est ainsi que nous allons procéder, en enregistrant la sélection de tous les enregistrements dont le code postal ne reflète pas une appartenance à la France outre mer.

Ouvrez la table d’attributs (clic droit sur la couche dans l’explorateur de couches). En bas à gauche de l’explorateur sélectionnez l’entrée déroulante Filtre avancé (Expression)

QGIS ouvre une fenêtre qui permet de construire une expression qui va permettre de filtrer les enregistrements.

Dans la liste des fonctions, sélectionnez le champ code_postal_uai dans Champs et valeurs, puis l’opérateur LIKE dans Opérateurs, puis tapez le format de valeur ‘97%’ qui signifie que l’on souhaite une valeur qui débute par les caractères ‘9’ et ‘7’ suivis de n’importes quels caractères. Vous devriez être familier du format d’expression si vous connaissez SQL, QGIS utilise un format semblable avec un jeu d’instructions SQL limité pour les expressions. Il suffit de chercher les enregistrements dont le code postal commence par 97 ou 98 et inverser la condition de recherche.

Le format complet de l’expression est :

NOT (
 "code_postal_uai"  LIKE '97%'
 OR
 "code_postal_uai"  LIKE '98%'
)
Une fois la table filtrée, sélectionnez tous les enregistrements par un clic dans l’angle gauche de la table

 

Il reste à enregistrer ces informations comme une nouvelle couche. Vous pouvez au passage vérifier sur la carte : les points sélectionnés apparaissent en jaune.

Faites un clic droit sur la couche d’établissements dans l’explorateur de couches tout en conservant la table d’attributs ouverte, sélectionnez l’entrée de menu Sauvegarder sous

QGIS vous propose de sauvegarder la couche, nous pourrions profiter de l’occasion pour transformer la couche dans une autre projection et basculer l’encodage des attributs à UTF-8 mais ne le faites pas pour le moment.

Sélectionnez le format ESRI Shapefile, le SCR de la couche RGF93 / Lambert-93, le codage ISO-8859-1 et cochez la case qui permet de n’enregistrer que les entités sélectionnées et celle qui permet d’ajouter le fichier sauvegardé à la carte. Validez puis supprimez la couche d’établissements obsolète.
Notez que QGIS limite la taille des noms d’attributs à 10 caractères lors de la sauvegarde : c’est une limitation du format Shapefile où les noms d’attributs sont limités à 10 caractères, les noms d’attribut du fichier des établissements sont donc tronqués.

 

Si vous parcourez les attributs de la couche vous pouvez constater que l’import du fichier texte ne s’est pas déroulé comme l’on pouvait s’attendre, le programme a considéré le champ code postal comme une valeur numérique et a supprimé le zéro qui préfixe les codes postaux sur quatre chiffres ! Nous allons rétablir le code postal en tant qu’attribut de type texte, complété à gauche par des caractères zéro ‘0’ sur 5 caractères.

Ouvrez les propriétés de la couche sauvegardée d’un double clic sur son nom, et sélectionnez l’onglet Champs.Cliquez sur l’icône de crayon pour basculer en mode édition : vous pouvez alors accéder à l’icône de boulier qui permet de calculer un champ.Ouvrez le calculateur de champ.

Le calculateur de champ permet de créer de nouveaux champs à partir de fonctions ou de champs existants, ou de redéfinir des champs existants. Ce qui nous importe est de redéfinir le champ de code postal afin qu’il devienne de type chaîne de caractère et soit complété à droite.

Cochez la case de création d’un nouveau champ nommé code_post de type texte et de longueur 254.
Dans la liste des fonctions vous pouvez trouver les fonctions de chaîne de caractères et plus bas les champs et valeurs, ces fonctions permettent à nouveau de construire une expression.Saisissez l’expression :
lpad("code_posta" , 5, '0')
L’expression signifie que l’on créé un nouveau champ en appliquant la fonction de chaîne de caractères complémentée à gauche sur le champ code_posta (code_postal_uai du fichier d’origine mais tronqué à 10 caractères) paramétrée pour complémenter avec des caractères ‘0’ sur une longueur de 5 caractères au total.Sélectionnez la colonne code_posta obsolète et cliquez sur l’icône de suppression à gauche du bouton de passage en édition puis désactivez le mode édition et enregistrez lorsque vous y êtes invité.

Vérifiez la modification dans la table d’attributs.

Nous allons maintenant procéder à une jointure spatiale entre deux couches de données afin de palier au manque d’informations de nos données d’établissement.

Effectuer une jointure entre 2 couches sous QGIS

Pour faire une jointure nous allons nous baser sur la position géographique des établissements et les données des fichiers IGN pour les contours administratifs des départements et communes.

Ajoutez la couche vecteur des départements récupérée lors de la préparation des données. (Faites glisser le fichier DEPARTEMENT.SHP sur la fenêtre d’application QGIS)
Rendez vous dans le menu Vecteur > Outils de gestion de données > Joindre les attributs par localisation
Paramétrez la fenêtre :

  • Indiquez une couche vecteur : indiquez la couche des établissements sur laquelle va se faire la jointure
  • Joindre la couche vecteur : indiquez la couche des données à joindre, les départements
  • Résumé de l’attribut : indiquez de prendre les attributs de la première entité au cas où plusieurs entités sont concernées par la jointure
  • Saisissez un fichier Shapefile de résultat
  • Table en sortie : lorsque les entités de la couche de départ n’ont pas de correspondance dans la couche à joindre on conserve tout de même les enregistrements

 

Après vérification dans la table attributaire, l’opération a aboutit à l’ajout des attributs de la couche département aux attributs de la couche d’établissements sauf pour une vingtaine d’établissements en bord de littoral. Cette fois recommençons l’opération en supprimant les entités de la table en sortie afin d’éliminer ces cas que nous négligerons donc parmi environ 65000 établissements.

QGIS utilise un opérateur intersection pour la jointure spatiale, je n’ai pas testé le cas où plusieurs entités seraient concernées par la jointure

Procéder à l’import

Le moment est venu d’importer les données d’établissements et départements en base. L’extension DB Manager (installée par défaut) va nous permettre de transférer ces couches vers PostGIS.

La première chose à faire avant de lancer le gestionnaire de bases de données est d’établir une connexion à la base PostGIS.

Ouvrez la fenêtre de gestion des tables PostGIS, depuis le menu
Couche > Ajouter une couche PostGIS

La fenêtre qui s’ouvre permet d’ajouter une table spatiale récupérée depuis PostGIS, c’est également à partir de cet endroit que l’on peut créer une nouvelle connexion, l’extension DB Manager ne permet pas de le faire contre toute attente.

Paramétrez une nouvelle connexion comme sur l’exemple suivant.Donnez un nom parlant à la connexion : je donne l’adresse du serveur, le type PostGIS et le nom de la base de données pour les distinguer plus tard depuis le nom.

QGIS viens d’enregistrer une nouvelle connexion à la base. Vous pouvez ensuite quitter la fenêtre d’ajout de tables PostGIS, nous reviendrons plus tard après création des tables.

Sous QGIS, ouvrez le gestionnaire de base de données :
Base de donnée > Gestionnaire de base de données > Gestionnaire de base de données

Une fois que vous avez la connexion vous pouvez charger les données d’une couche ou d’un fichier à partir du gestionnaire de base de données.

Sélectionnez la connexion dans l’explorateur, sous l’entrée PostGIS, puis cliquez sur l’icône en forme de flèche vers la gauche, elle permet de paramétrer un import.

La couche des départements est dans le système de coordonnées de référence EPSG:2154 – RGF93 / Lambert-93, l’encodage est System. (Ces informations sont visibles dans les propriétés générales de la couche). Pour importer les données nous préciserons de convertir l’encodage en UTF-8 et de passer le SCR des données en WGS 84 d’identifiant EPSG:4326.

Ce dernier point n’est pas obligatoire mais je préfère unifier les système pour faciliter les développements ultérieurs et éviter les conversions en aval entre tables de données sous PostGIS.

Paramétrez l’import de la couche des départements :

  • Saisie : vous pouvez sélectionner une couche parmi les couches du projet ouvert ou bien utiliser le bouton de navigation à droite afin de préciser le chemin vers un fichier Shapefile.
    Sélectionnez la couche DEPARTEMENT
  • Saisissez le nom de la table Departement pour la table en sortie
  • Précisez le schéma : public
  • Dans l’encart Action vous pouvez si besoin écraser la table existante
  • Précisez la clef primaire : id_geofla. Ici on réutilise l’attribut qui sert d’identifiant dans la couche, par défaut le gestionnaire créé un identifiant associé à une séquence si ce champ n’est pas renseigné
  • Précisez si vous le souhaitez la colonne de géométrie the_geom. Certaines conventions utilisent the_geom, depuis les dernières versions de PostGIS geom est utilisé comme nom de colonne par défaut
  • Précisez le SCR source : 2154
  • Précisez le SCR cible : 4326
  • Précisez le codage cible : UTF-8
  • Cochez la case afin de créer un index spatial

Une fois l’import réussi il est possible de réactualiser la liste des objets de la connexion afin de faire apparaître la table qui vient d’être crée. Les onglets du gestionnaire sur la droite permettent d’accéder aux informations sur la table, de vérifier le contenu et le rendu.

Le gestionnaire donne également accès à certaines vues, si l’on consulte la vue geometry_columns on vérifie bien la présence d’une nouvelle ligne où l’on accède aux informations sur la colonne de type Geometry que l’on a créé dans la nouvelle table. Vous pouvez à nouveau vous référer aux documents de standards sur la vue geometry_columns. La vue référence bien la colonne geom de la table Departement, de type MULTIPOLYGON et de srid 4326.

Ultérieurement j’ai rencontré des problèmes avec des requêtes spatiales sur cette table du fait que la colonne est de type MULTIPOLYGON, mais que certaines de ses entrées sont de type GeometryCollection. Une solution serait peut être de modifier dans les fichiers SQL générés les instructions de création de table. J’ai choisi par manque de temps d’ignorer (grâce à la fonction ST_IsValid()) les lignes de type GeometryCollection qui représentent un pourcentage négligeable.

Si vous observez maintenant la table sous PgAdmin vous pouvez vérifier le contenu de la colonne de géométrie sur quelques enregistrements : le contenu est au format standard WKB (Well-Known Binary), un format de représentation hexadécimal. A présent vous pouvez parcourir les fonctions spatiales de PostGIS dans l’explorateur d’objet de PgAdmin, la fonction st_astext(geometry) permet de retourner la représentation de la colonne dans le format standard WKT (Well-Known Text) que l’on utilisera ultérieurement pour les échanges d’informations entre les différents composants de l’application cartographique que l’on va mettre en place.

Un exemple de requête retournera une représentation WKT où l’on va voir apparaître des objets de type MULTIPOLYGON, le type de la colonne, avec des coordonnées en 2 dimensions (X et Y) exprimées dans le SCR de référence SRID 4326 tel que l’on peut le consulter dans la définition de la colonne dans la vue geometry_columns :

SELECT st_astext(geom) FROM "Departement" LIMIT 1;
"MULTIPOLYGON(((5.83122641362104 45.9384595782932,5.82212102507634 45.9301359535238,5.82915456688181 45.9172371235251,5.82614026745121 45.903273829068,5.81515407306268 45.8772007832312,5.80752504239684 45.8606398467483,5.80110949497983 45.8457964776072,5.79 (...)"

Procédez à l’import de la couche des établissements via DB Manager, nommez la table Etablissement et utilisez l’attribut numero_uai comme clef primaire.(Mon installation sur VPS à ressources limitées ne supporte pas le chargement de telles quantités de données au travers de la connexion au serveur de base de données. Pour les développements je travaillerai en local pour le confort d’utilisation et j’utiliserai l’utilitaire de sauvegarde et restauration pour migrer la base afin de faire fonctionner le prototype)

Attention à l’encodage des caractères… malgré les options dans DB Manager j’ai dû au préalable changer le jeu de caractères en UTF-8 en sauvegardant la couche sous un autre nom (on peut préciser le jeu de caractères lors de l’enregistrement)

Pour gagner du temps je n’ai pas supprimé d’attributs inutiles ni renommé les noms tronqués des colonnes… mais rien ne vous empêche de le faire. Lors de la sauvegarde d’une couche on peut renommer les attributs, sinon après import il suffit de passer par PgAdmin et d’effectuer un clic droit sur la table dans l’explorateur d’objets pour faire afficher les propriétés. On accède alors à la définition des colonnes dans un onglet, où l’on peut ajouter, modifier ou supprimer une colonne.

Import des données de parcelles

Traiter les parcelles une à une serait assez consommateur de temps.

Nous disposons d’une centaine de ressources sous format d’archive zip, d’en moyenne 15M, ce qui totalise environ 1,5Go.

Je fais le choix de conserver le découpage par départements plutôt que d’utiliser une table unique où l’on ajouterai les données de chaque département. Cela permettra d’interroger individuellement les tables en fonction des codes département.

Deux options se présentent alors pour interroger la France métropolitaine dans son intégralité : dans un premier temps j’utiliserai une vue, puis je couvrirai les capacités de partitionnement de tables offertes par PostgreSQL avec héritage d’une table maître pour comparer les approches.

Procéder à l’import

L’import doit se dérouler de manière automatique, on doit pour cela utiliser un script combiné à un utilitaire d’import en ligne de commande.

Pour le langage de script j’utiliserai Python qui est notamment employé sous QGIS.

Vous pouvez trouver une version embarquée sous Windows sous le répertoire d’installation de QGIS :

cd C:\Program Files\QGIS Chugiak\bin
python --version
Python 2.7.4

Pour l’utilitaire en ligne de commande j’utiliserai shp2pgsql pour la génération de fichiers SQL couplé à psql pour le chargement en base, il sont disponibles dans le répertoire d’installation de PostgreSQL :
C:\Program Files\PostgreSQL\9.3\bin

Ajouter le chemin vers ces exécutables à la variable système PATH permet de s’affranchir des chemins lorsque l’on lance une commande.

Import à l’aide de shp2pgsql

shp2pgsql convertit un fichier Shapefile en instructions SQL qui peuvent être exploitées ensuite pour alimenter une base de donnée, par psql dans notre exemple.

L’utilitaire peut être accompagné d’une interface utilisateur qui permet également le chargement (ou l’export) de plusieurs fichiers. Elle offre cependant moins de contrôle et je ne couvrirai pas son utilisation mais elle peut s’avérer pratique lorsque l’on a un nombre restreint de fichiers à importer.

Sous Windows l’utilitaire est disponible sous le répertoire d’installation de PostgreSQL :

C:\Program Files\PostgreSQL\9.3\bin\postgisgui\shp2pgsql-gui.exe

La sortie de commande de shp2pgsql peut être capturée dans un fichier SQL (ou redirigée vers un autre utilitaire sous système de type Unix)

Les paramètres de commande sont :

shp2pgsql [<options>] <shapefile> [[<schema>.]<table>]

L’ouverture d’un fichier *.shp sous QGIS nous renseigne sur le SCR et l’encodage. Un coup d’œil aux attributs permet d’ignorer les problématiques d’encodage puisque l’on travaille avec des valeurs numériques.

Au niveau des options dans notre cas de figure :

  • -s [:] Fixe le SRID en entrée et en sortie
  • -c Créé une nouvelle table et les instructions d’insertion
  • -I Créé un index spatial sur la colonne géométrie
  • -N skip Les enregistrements avec géométries vides ne seront pas importés
Si vous ne l’avez pas déjà fait, décompressez les archives de parcelles dans un répertoire !

L’usage est alors le suivant :

(Sous Windows 7 la combinaison touche SHIFT + Clic droit sur la fenêtre d’explorateur en cours permet d’ouvrir une fenêtre de commande à cet endroit, vous pouvez sinon ouvrir un terminal à l’aide de l’utilitaire cmd)

cd C:\SIG\data\parcelles\RPG_2012_004
shp2pgsql -s 2154:4326 -c -I -N skip RPG_2012_004.shp Cultures_004 > "C:\SIG\sql\Cultures_004.sql"

Le fichier SQL généré contient dans une transaction les instructions de création de table, d’indexe et de clef primaire, ainsi que les instructions d’insertion des données. La colonne géométrie est ajoutée après création de table pour être compatible avec les anciennes versions de PostGIS. Nous pourrons plus tard altérer ces instructions SQL pour tester l’implémentation avec héritage d’une table maître.

SET CLIENT_ENCODING TO UTF8;
SET STANDARD_CONFORMING_STRINGS TO ON;
BEGIN;
CREATE TABLE "cultures_004" (gid serial,
"num_ilot" varchar(12),
"cult_maj" int4);
ALTER TABLE "cultures_004" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('','cultures_004','geom','4326','MULTIPOLYGON',2);
INSERT INTO "cultures_004" ("num_ilot","cult_maj",geom) VALUES ('004-190399','4',ST_Transform('01060000206A0800000100000001030000000100000015000000787AA5EC64362C41C9E53F4C6D195841A8A44E8015362C41BC7493087E1958417DAEB622F9352C4106F016707D195841151DC92599352C410DE02D509E195841F241CFA696352C41F697DD03A01958416519E2189C352C41D2DEE0DBA119584183C0CAE1B1352C411A2FDDE4A4195841E10B9369BB352C41EA73B565A6195841151DC965C5352C41840D4FB7A9195841AACFD556CD352C41F1F44A19AC1958416D348037E2352C4189B0E139AE19584182E2C71812362C4147E17AE4B0195841C56D34405F362C41627FD9EDA7195841499D80E673362C415E29CB20A4195841B515FB4B7C362C41D634EF20A01958416132557078362C41779CA2338E19584179C729FA7C362C41D1915C7E85195841174850BC86362C41E92631F87E195841F163CCDDA7362C41F54A590E73195841728A8EA487362C41627FD99D70195841787AA5EC64362C41C9E53F4C6D195841'::geometry, 4326));
...
CREATE INDEX "cultures_004_geom_gist" ON "cultures_004" USING GIST ("geom");
COMMIT;

Il reste à exécuter le script sur la base de données. Nous pouvons utiliser l’utilitaire psql.

psql --host=localhost --port=5432 --username=postgres --no-password --dbname=Pesticides --file="C:\SIG\sql\Cultures_004.sql"

Les options sont :

  • –host=HOTE nom d’hôte du serveur de la base de données ou répertoire de la socket (par défaut : socket locale)
  • –port=PORT port du serveur de la base de données (par défaut : « 5432 »)
  • –username=NOM nom d’utilisateur de la base de données
  • –no-password ne demande jamais un mot de passe
  • –dbname=NOM_BASE indique le nom de la base de données à laquelle se connecter
  • –file=FICHIER exécute les commandes du fichier, puis quitte

Vous pouvez vérifier l’import sous PgAdmin et supprimer la table en cascade.

Nous savons comment réaliser un import, voyons comment réaliser un traitement par lot.

Import automatisé, traitement par lot à partir d’un script

Un script très simple va permettre d’automatiser l’appel à l’utilitaire shp2pgsql avec les paramètres que nous avons validés précédemment.

Ajoutez si ce n’est pas déjà fait le chemin vers l’exécutable python à la variable système PATH

Enregistrez le script suivant dans un fichier d’extension .py placé sous le répertoire où se trouvent vos données, le script utilise le répertoire courant pour rechercher les fichiers Shapefile.

J’ai pris soin d’écrire du code portable mais vous devrez peut être procéder à des adaptations sous votre système d’exploitation.

#! /usr/bin/python
import os
import fnmatch

instructions = r'shp2pgsql -s 2154:4326 -c -I -N skip "%s" %s > "%s"'

currentdir = os.path.realpath(os.curdir)
sqldirectory = currentdir + os.sep + 'sql'
if not os.path.exists(sqldirectory):
    os.makedirs(sqldirectory)

for root, dirnames, filenames in os.walk(currentdir):
  for shapefilename in fnmatch.filter(filenames, '*.shp'):
      print shapefilename, '...'
      shapefilepath = root + os.sep + shapefilename
      tablename = 'Cultures_' + shapefilename[-6:-4]
      sqlfilepath = sqldirectory + os.sep + tablename + '.sql'
      command = instructions % (shapefilepath, tablename, sqlfilepath)
      os.system(command)

raw_input("Appuyez sur une la touche Entree pour quitter...")

Que fait le programme ?

Premièrement il créé une variable où l’on conserve un gabarit d’appel générique à la commande shp2pgsql, le nom du fichier shapefile, le nom de table et le fichier de destination sont des paramètres.

Ensuite il récupère le répertoire courant et créé si besoin un répertoire sql sous le répertoire courant.

Le programme parcours alors le répertoire courant et les sous répertoires à la recherche de fichiers d’extension .shp, lorsqu’il trouve un fichier il extrait le code du département et construit les paramètres (chemin vers le fichier, nom de table, chemin vers le fichier de destination) et appelle la commande ainsi formatée en exécution système.

Nos fichiers sql sont alors disponibles sous le répertoire sql dans le répertoire courant, il reste à les modifier à nos besoins et à les importer.

Ce deuxième script est très semblable au premier. Dans ce cas on recherche les fichiers d’extension .sql sous les dossiers du répertoire courant et lorsqu’ils sont trouvés la commande système d’appel à psql est lancée afin d’exécuter les instructions SQL contenues dans le fichier.

#! /usr/bin/python
import os
import fnmatch

instructions = r'psql --host=localhost --port=5432 --username=postgres --no-password --dbname=Pesticides --file="%s"'

for root, dirnames, filenames in os.walk(os.path.realpath(os.curdir)):
  for sqlfilename in fnmatch.filter(filenames, '*.sql'):
      print sqlfilename, '...'
      sqlfilepath = root + os.sep + sqlfilename
      command = instructions % (sqlfilepath)
      os.system(command)

raw_input("Appuyez sur une la touche Entree pour quitter...")

L’opération prend un certain temps, même en local, mais le chargement en base c’est bien déroulé.

A présent créons une vue qui regroupera l’information de toutes les tables, ou plutôt une vue matérialisée.

Une vue matérialisée, à la différence d’une vue, stocke à la fois la requête qui a permis la génération de la vue mais également les données du résultat de requête. La vue matérialisée est équivalente sur ce point à une table et il est possible de l’indexer de la même façon. Ce comportement diffère de la vue qui exécute la requête à chaque appel, ce qui dans notre cas se traduirait par des performances moindres, dans le sens où les indexes de chaque table ne seraient probablement pas utilisés au mieux.

L’inconvénient principal des vues matérialisées c’est qu’il faut rafraîchir la vue lorsque les tables impactées par la requêtes sont modifiées, ce qui ne pose pas problème par rapport à nos données qui sont des données destinées à la consultation uniquement, pas à la mise à jour. Par contre les données sont dupliquées, nous aurons donc environ 1,5Go de données en double.

Lors de la création des instruction SQL, shp2pgsql créé une colonne gid qui sert de clef primaire. Cette colonne est construite sur la base d’une séquence à partir de 1, incrémentée de 1, ce qui est problématique pour lier plusieurs tables dans une vue : nous ne devons pas utiliser l’attribut gid dans la vue, ce qui créerait des doublons de clef, mais préciser lors du chargement d’une couche sous QGIS que l’attribut de référence pour la clef est le numéro d’îlot num_ilot

Nous devons aboutir à une instruction SQL comme celle ci :

CREATE MATERIALIZED VIEW mv_cultures AS
SELECT num_ilot, cult_maj, geom FROM "cultures_01" UNION ALL
SELECT num_ilot, cult_maj, geom FROM "cultures_02" UNION ALL
(...)
SELECT num_ilot, cult_maj, geom FROM "cultures_95"

Nous possédons tous les codes département dans la table Departement, pour générer la vue il va falloir créer la requête par une instruction SQL qui récupère tous les codes département sauf les Hauts de seine et Paris pour lesquels il n’existe pas de table, puis exécuter la requête :

SELECT
 'CREATE MATERIALIZED VIEW mv_cultures AS ' ||
 string_agg(format('SELECT num_ilot, cult_maj, geom FROM "cultures_%s"', lower("code_dept")), ' UNION ALL ')
FROM
 "Departement"
WHERE
 code_dept NOT IN ('75', '92');

Vous pouvez enregistrer la requête dans un fichier et passer par psql pour générer le résultat et vérifier avant exécution.

psql --host=localhost --port=5432 --username=postgres --no-password --dbname=Pesticides --file="C:\SIG\sql\createquery.sql" --quiet --output="C:\SIG\sql\query.sql" --log-file="C:\SIG\sql\createquery.log"
psql --host=localhost --port=5432 --username=postgres --no-password --dbname=Pesticides --file="C:\SIG\sql\query.sql"

Procédez de même avec une instruction CREATE VIEW et une instruction CREATE TABLE afin de comparer les plans d’exécutions de requêtes de type SELECT.N’oubliez pas de créer clefs primaires et indexes sur la vue matérialisée et la table…

CREATE INDEX t_cultures_geom_gist ON t_cultures USING gist (geom);
ALTER TABLE t_cultures ADD PRIMARY KEY (num_ilot);

CREATE INDEX mv_cultures_geom_gist ON mv_cultures USING gist (geom);
CREATE UNIQUE INDEX ON mv_cultures (num_ilot);
Partitionnement de la table

Une autre option qui s’offre à nous est de créer une table logique à partir de plusieurs tables physiques plus petites.

PostgreSQL supporte une forme simple de partitionnement de table. Je vais l’essayer sur ce cas de figure pour le découpage départemental des données et comparer le plan d’exécution avec les autres approches.

Le fichier de configuration de PostgreSQL permet de contrôler de quelle façon le planificateur de requêtes utilise les contraintes de table pour optimiser les requêtes.

Pour modifier le paramètre d’exclusions de contraintes, il suffit d’éditer le fichier postgresql.conf (sous C:\Program Files\PostgreSQL\9.3\data) :

constraint_exclusion = partition	# on, off, or partition

Le paramètre par défaut est partition, il pourra être intéressant de modifier le paramètre pour vérifier l’impact sur les plans d’exécution.

Procédons à la création de la table logique qui sera la table maître dont les autres tables vont hériter. La table ne possède pas d’indexe, et par rapport à la définition des tables cultures, elle possède désormais un code département.

CREATE TABLE master_cultures
(
  num_ilot character varying(12),
  cult_maj integer,
  code_dept character varying(2),
  geom geometry(MultiPolygon,4326)
)
;

Il reste ensuite à créer chaque table fille avec les instructions d’héritage appropriées. Les requêtes sont construites à la volée comme précédemment, le code n’est pas très élégant mais sera à usage unique !

SELECT
 string_agg(
	 format('DROP TABLE IF EXISTS child_cultures_%s;CREATE TABLE child_cultures_%s AS SELECT num_ilot, cult_maj, (SELECT CODE_DEPT from "Departement" WHERE code_dept = ''%s''), geom FROM "cultures_%s";ALTER TABLE child_cultures_%s ADD CONSTRAINT CK_%s CHECK (code_dept = ''%s'');ALTER TABLE child_cultures_%s ADD CONSTRAINT PK_%s PRIMARY KEY (num_ilot);CREATE INDEX child_cultures_%s_geom_gist ON child_cultures_%s USING gist (geom);ALTER TABLE child_cultures_%s INHERIT master_cultures;'
	 , lower("code_dept"), lower("code_dept"), lower("code_dept"), lower("code_dept"), lower("code_dept"), lower("code_dept")
	 , lower("code_dept"), lower("code_dept"), lower("code_dept"), lower("code_dept"), lower("code_dept"), lower("code_dept")
 ), ' ; ')
FROM
 "Departement"
WHERE
  code_dept NOT IN ('75', '92');

Ces données ne sont pas destinées à être mises à jour, il n’est pas nécessaire de contrôler les opérations de mise à jour sur la table maître par des procédures stockées comme dans la documentation. (La table ne possède d’ailleurs pas de clef primaire ni d’OIDS, l’accès en édition sera interdit sous PgAdmin)

Tout ça pour ça… Et maintenant ? Contrôlons la validité de chaque vue et chaque table créées pour un usage en tant que couche vectorielle.

Contrôle du chargement des couches PostGIS sous QGIS

Afin de vérifier le chargement d’une couche PostGIS à partir de la vue, de la vue matérialisée ou de la table il faut au préalable définir une emprise pour le rendu qui soit suffisamment restreinte pour ne pas charger toutes les données depuis PostGIS sans quoi il va falloir être à nouveau très patients !

Depuis la couche GEO_FLA des communes j’ai effectué sous QGIS un zoom sur Montélimar.

Cela se fait simplement depuis la table d’attributs en sélectionnant une ligne puis en cliquant sur l’icône de zoom sur la ligne sélectionnée.

La barre de statuts en bas permet d’afficher l’emprise (icône souris à côté des coordonnées), elle s’exprime sous la forme xmin,ymin : xmax,ymax où les coordonnées sont exprimées dans le SCR du projet :

xmin,ymin : xmax,ymax EPSG:4326
4.689,44.509 : 4.803,44.601
Sous QGIS, zoomez sur une commune pour restreindre les limites d’affichage puis allez dans le menu :
Couche > Ajouter une couche PostGISSous la connexion à la base vous devez apercevoir les vues et table nouvellement crées.

Avant d’ajouter une vue, vous devez préciser à QGIS la clef primaire, le numéro d’îlot, sans quoi QGIS retournera une erreur d’invalidité de la couche.

L’ajout de la couche à partir de la table ou de la vue matérialisée est presque instantané, alors que pour la vue… eh bien, vous vous féliciterez peut être d’avoir investi dans un disque SSD dernier cri.

A l’œil nu, le rafraîchissement des polygones sur la carte prend sensiblement le même temps pour la table et la vue matérialisée lorsque l’on agrandit l’étendue représentée, il est à nouveau beaucoup plus lent en ce qui concerne la vue.

A partir de la table maître et des tables héritées l’ajout d’une couche n’est pas instantané, une douzaine de secondes sont nécessaires pour initialiser le rendu par contre le rafraîchissement est ensuite beaucoup plus réactif que pour la table et la vue matérialisée lorsque l’on représente une plus grande étendue.

Cherchons à savoir pourquoi les performances sont mauvaises avec une simple vue, quelles optimisations apporter aux autres cas et quel choix faire pour notre application.

Bilan à l’heure du choix. Explication du plan d’exécution

PostgreSQL permet de visualiser les plans d’exécution de chaque requête SQL, l’instruction EXPLAIN retourne des informations pertinentes sur le plan d’exécution qu’il va utiliser.

Pour tester les plans d’exécution dans chaque configuration je vais construire une requête spatiale afin de simplement retourner toutes les parcelles sur une étendue géographique donnée.

En effet, ce sera le mode de fonctionnement de l’application web, il importe de le tester : les librairies javascript en charge du rendu des cartes vont faire appel aux informations sur l’étendue de la carte en cours, avec possibilité d’agrandir ou rétrécir l’étendue représentée.

Lorsque l’étendue sera élevée la couche des parcelles ne sera simplement pas affichée, car elle n’apporte pas d’information particulière dans un contexte de rendu global, par contre elle sera visible dès que l’échelle de la carte sera suffisamment élevée.

Sous PostGIS, la fonction ST_MakeEnvelope créé un polygone rectangulaire à partir de coordonnées X, Y et d’un SRID.

Le format est le suivant :

geometry ST_MakeEnvelope(double precision xmin, double precision ymin, double precision xmax, double precision ymax, integer srid=unknown);

Pour la commune de CAHUZAC située dans le département de code 47, le polygone recouvrant l’étendue géographique sera donc obtenu par :

SELECT ST_AsText(ST_MakeEnvelope(0.533, 44.639, 0.589, 44.685, 4326));

L’option ANALYSE permettra de lancer l’exécution de la requête afin de disposer des temps d’exécution et nombre de lignes réels.

L’étendue géographique sera celle de l’Auvergne :

EXPLAIN ANALYSE SELECT ST_AsText(geom) FROM t_cultures WHERE geom && ST_MakeEnvelope(2, 45 , 5, 47, 4326);
EXPLAIN ANALYSE SELECT ST_AsText(geom) FROM mv_cultures WHERE geom && ST_MakeEnvelope(2, 45 , 5, 47, 4326);
EXPLAIN ANALYSE SELECT ST_AsText(geom) FROM v_cultures WHERE geom && ST_MakeEnvelope(2, 45 , 5, 47, 4326);
EXPLAIN ANALYSE SELECT ST_AsText(geom) FROM master_cultures WHERE geom && ST_MakeEnvelope(2, 45 , 5, 47, 4326);

Les résultats sur cette requête me laissent perplexe car ils semblent contradictoires avec le ressenti sous QGIS.

La méthode est certainement sujette à un problème, si quelqu’un peut m’éclairer les remarques seront bienvenues pour construire des requêtes de test plus pertinentes !

Ce article est assez long, et je m’égare. La prochaine étape sera d’installer un serveur cartographique. L’héritage de table donne des résultats corrects même si la vue semble meilleure, par rapport au ressenti sous QGIS c’est la solution que je vais retenir par défaut, je conserve toutefois la possibilité de basculer sur les autres solutions.

Installation d’une base de données spatiales PostgreSQL PostGIS

PostgreSQL PostGIS

Introduction

PostgreSQL est un système de gestion de bases de données relationnelles Open Source. Il est reconnu pour sa fiabilité et ses performances. PostGIS est un module de PostgreSQL qui apporte un lot de capacités et de fonctions spatiales à la base. Il respecte les standards internationaux et il est devenu une référence pour stocker, gérer et analyser des données spatiales. PostGIS est supporté par la fondation OSGeo.

PostGIS implémente les spécifications SQL pour l’information géographique OpenGIS® définies par l’Open Geospatial Consortium (OGC), un consortium international qui établie les documents techniques standards OGC® Standards. Ce document technique est une source d’informations qui peut servir même sans rentrer dans une étude complète et ennuyeuse du standard !

Vous pouvez suivre le tutoriel d’introduction à PostGIS proposé en anglais par OpenGeo. La documentation PostGIS propose des exemples standards des fonctions SQL, gardez là également sous la main !

Installation sous Debian « Wheezy »

Une ligne de commande suffit sous Debian pour installer le serveur de base de données et l’extension spatiale cependant nous allons au préalable ajouter un dépôt au gestionnaire de paquets de la distribution Debian Wheezy afin de récupérer les dernières versions de Postgresql et PostGIS.

Éditez le fichier de sources du gestionnaire de paquets
vi /etc/apt/sources.list.d/pgdg.list
Insérez la ligne suivante :
deb http://apt.postgresql.org/pub/repos/apt/ wheezy-pgdg main
Récupérez la clef de sécurité et mettez à jour les informations dans le gestionnaire de paquets
wget https://www.postgresql.org/media/keys/ACCC4CF8.asc
apt-key add ACCC4CF8.asc
apt-get update
Procédez à l’installation
apt-get install postgresql postgresql-contrib libgdal-dev postgresql-9.3-postgis-2.1

Vous pouvez vérifier le numéro de version précisément installée.

psql --version
psql (PostgreSQL) 9.3.5

Lors de l’installation le système a créé un utilisateur système postgres, le client postgresql psql que nous venons d’utiliser est disponible dans le chemin des exécutables, et un super utilisateur nommé postgres, l’utilisateur par défaut des bases de données a été créé dans PostgreSQL.

Nous allons modifier le mot de passe de l’utilisateur base de données postgres. Notez le mot de passe, nous en aurons besoin ultérieurement.

su - postgres
psql
postgres=# \password postgres
Saisissez le nouveau mot de passe :
Saisissez-le à  nouveau :
postgres=# \q

Il s’agit de procéder rapidement à l’installation d’un serveur privé suffisant pour développer une application, je vais ignorer les aspects sécurité et performance d’un environnement de production, c’est à dire ne pas créer de rôles et d’utilisateurs de base de données et ne pas paramétrer le serveur au delà du strict minimum nécessaire à mon propre accès.

Dans l’immédiat ce qui m’importe est de pouvoir me connecter localement depuis l’outil en ligne de commande et depuis une application d’administration client distante avec, pour sécuriser l’accès, une restriction sur mon adresse IP et une connexion par mot de passe crypté.

PostgreSQL est paramétré par défaut pour refuser les connexions entrantes à partir des machines externes, ce qui implique que l’on ne peut pas se connecter depuis un outil d’administration sur un poste de travail.

Pour rendre l’accès disponible il faut paramétrer PostgreSQL pour écouter les connexions entrantes.

Paramétrez un accès distant :
su - postgres
vi /etc/postgresql/9.3/main/postgresql.conf
Modifiez le fichier pour écouter les connexions sur la boucle locale et sur l’adresse IP publique du serveur où est hébergé PostgreSQL.
#------------------------------------------------------------------------------
# CONNECTIONS AND AUTHENTICATION
#------------------------------------------------------------------------------

# - Connection Settings -

listen_addresses = 'localhost, 92.222.27.205'           # what IP address(es) to listen on;
                                        # comma-separated list of addresses;
                                        # defaults to 'localhost', '*' = all
                                        # (change requires restart)
port = 5432

Vous pouvez vérifier la prise en compte de vos modifications en lançant la commande netstat avant et après voir redémarré le serveur

exit
service postgresql restart
netstat -an | grep 5432
tcp        0      0 92.222.27.205:5432      0.0.0.0:*               LISTEN
tcp        0      0 127.0.0.1:5432          0.0.0.0:*               LISTEN
tcp6       0      0 ::1:5432                :::*                    LISTEN
...

Maintenant que PostgreSQL écoute les connexions depuis l’extérieur, il est nécessaire d’étendre les règles de contrôle d’accès par défaut pour accepter les connexions.

Ajoutez une règle d’accès au fichier pg_hba.conf qui autorisera de vous connecter avec mot de passe depuis l’adresse IP de votre poste de travail :
su - postgres
vi /etc/postgresql/9.3/main/pg_hba.conf
Ajoutez la ligne suivante, l’information host permet d’indiquer que la tentative de connexion doit être de type TCP/IP, la première information all spécifie que toutes les bases de données sont acceptées, la seconde information all permet d’accepter tous les noms d’utilisateur de la base de données, enfin la plage d’adresses permet de spécifier les machines client acceptées, saisissez votre adresse IP (l’hôte), suivie du masque CIDR /32 (255.255.255.255) afin d’accepter une seule machine hôte, le dernier paramètre md5 indique la méthode d’authentification, vous devrez vous connecter par mot de passe crypté md5 :
host all all 85.171.50.88/32 md5

Après avoir redémarré le serveur nous allons voir comment nous connecter depuis une application cliente installée sur notre poste de travail.

exit
service postgresql restart

Outil d’administration Pg Admin III

Pg Admin est un outil d’administration et de gestion pour PostgreSQL. C’est un logiciel libre disponible sous Windows, Mac, Linux. Nous utiliserons Pg Admin en priorité par rapport à l’outil en ligne de commande.

Enregistrer une nouvelle connexion

Téléchargez et installez Pg Admin III, sélectionnez dans l’entrée du menu principal
Fichier > Ajoutez un serveur

Une fenêtre s’ouvre ou vous devez paramétrer la connexion, saisissez le nom sous lequel sera enregistré la connexion, le nom d’hôte du serveur sur lequel tourne PostgreSQL, le nom d’utilisateur postgres et le mot de passe que vous avez renseigné après installation du serveur.

Félicitations, vous pouvez désormais vous connecter et utiliser le navigateur d’objet à gauche pour obtenir des informations et interagir avec PostgreSQL.

Créer une base de données et activer l’extension spatiale

Nous allons passer par Pg Admin pour la plupart des opérations d’administration de bases de données. Ces opérations pourraient très bien être réalisées en SQL sous l’outil de ligne de commande psql : Pg Admin vous donne accès au code SQL qu’il génère.

Dans l’explorateur d’objets, faites un clic droit sur l’objet Bases de données et sélectionnez l’entrée Ajouter une base de donnée dans le menu contextuel qui vient de s’ouvrir.Remplissez le formulaire d’ajout de base de données :

  • Dans l’onglet Propriétés donnez pour nom Pesticides, et pour propriétaire postgres
  • Dans l’onglet Définition sélectionnez le codage UTF8, la collation et le type de caractères fr_FR.UTF-8
  • Dans l’onglet SQL vous pouvez vérifier l’instruction SQL qui sera exécutée et au besoin décocher le mode lecture seule pour passer en édition

L’outil exécute la commande SQL sur la connexion serveur active :

CREATE DATABASE "Pesticides"
  WITH ENCODING='UTF8'
       OWNER=postgres
       LC_COLLATE='fr_FR.UTF-8'
       LC_CTYPE='fr_FR.UTF-8'
       CONNECTION LIMIT=-1;
Vous pouvez altérer certaines propriétés de la base ultérieurement d’un double clic sur l’écran de propriétés qui apparaît à la sélection de la base de données dans le navigateur d’objet mais prenez soin de bien définir à la création les propriétés relatives au jeu de caractères et à la localisation

Vous venez de créer une base de données, il vous reste à activer le module spatial pour cette base.

Dans l’explorateur d’objets, faites un clic droit sur l’objet base de données Pesticides et sélectionnez l’entrée Ajouter un objet > Ajouter une extension dans le menu contextuel qui vient de s’ouvrir.Sélectionnez le nom postgis

L’outil exécute la commande SQL d’ajout d’une extension à la base de données :

CREATE EXTENSION postgis;

L’activation de l’extension permet désormais d’utiliser les fonctions spatiales.

Remarquez la création d’une table spatial_ref_sys sous le schéma public (par défaut si vous n’avez pas modifié le schéma). Vous pouvez parcourir les 100 premières lignes à l’aide d’un clic droit sur la table dans l’explorateur d’objet, sélectionnez l’entrée de menu contextuel Afficher les données > Visualiser les 100 premières lignes.

Vous pouvez vérifier dans la documentation du standard OpenGIS® les informations relatives à cette table : elle sert à décrire le système de coordonnées et les transformations pour la représentation de la géométrie. Elle permet de conserver en base les informations sur chaque système spatial de référence :

  • SRID est l’identifiant du système spatial de référence, il est donc unique et sert de clef dans la table
  • AUTH_NAME est le nom de l’autorité qui enregistre le système spatial, typiquement EPSG (European Petroleum Survey Group) pour le jeu de données par défaut. (voir http://spatialreference.org/ pour d’autres systèmes)
  • AUTH_SRID est l’identifiant du système de référence au sein de l’autorité qui l’enregistre (le code de projection EPSG dans le cas d’EPSG)
  • SRTEXT est la description du système spatial de référence au format Well known text (WKT est un format texte standard pour représenter les objets géométriques et les informations rattachées)
  • La dernière colonne PROJ4TEXT n’est pas dans le standard. PostGIS utilise la librairie de projections cartographiques Proj4 pour les transformations de coordonnées, la colonne contient les définitions de coordonnées Proj4 pour le SRID
A retenir : La table spatial_ref_sys contient les informations nécessaires pour permettre d’identifier les systèmes spatiaux de référence (SRS) et permettre de transformer et changer de projection d’un système à un autre.

Créer une table spatiale se résume à ajouter une colonne de type Geometry à la définition de table où l’on va préciser certaines informations en plus des informations de colonne par défaut :

  • Le SRID du système spatial de référence qui servira de clef étrangère vers la table spatial_ref_sys
  • La dimension spatiale COORD_DIMENSION : indique un entier qui correspond à 2, 3 ou 4 dimensions pour les coordonnées
  • Le type de l’objet géométrique qui doit être enregistré si l’objet possède un type unique (POINT, LINESTRING, POLYGON, MULTIPOINT, MULTILINESTRING, MULTIPOLYGON, GEOMETRYCOLLECTION, …) ou bien le type GEOMETRY

Nous n’allons pas créer manuellement les tables spatiales mais munis de ces informations nous pouvons passer au chargement de données et contrôler le format des tables crées.

Préparation des données sous QGIS Desktop

Récupération des données nécessaires à la mise en oeuvre de l’application cartographique

Vous pouvez récupérer la plupart des données que nous allons exploiter depuis de la plateforme ouverte des données publiques françaises sur data.gouv.fr

Données et coordonnées géographiques des établissements scolaires français

Où se procurer les données ?

L’Education Nationale a mis à disposition sur la plateforme data.gouv.fr un fichier d’adresses et géolocalisation des établissements d’enseignement du premier et second degrés. (Le premier degré correspond aux écoles maternelles et élémentaires, le second degré regroupe collèges et lycées)

Téléchargez la version txt des données d’adresses, notez que le fichier est récent, il date de juin 2014. Prenez également connaissance de la notice explicative pour les données d’adresses.

Le fichier possède de nombreuses informations exploitables, premièrement le numéro de l’unité administrative immatriculée (uai). Ce numéro permet d’identifier un établissement, il servira d’identifiant ultérieurement lorsque nous chargerons les établissement dans une base de données. Au travers de ce numéro nous pourrons également rediriger l’utilisateur vers la base centrale des établissement de l’éducation nationale où il pourra trouver les informations sur un établissement dont nous ne disposons pas dans le fichier.

La base centrale des établissement permet d’effectuer une recherche par numéro uai, mais faute de données sur le nombre d’élèves par établissement l’application ne pourra pas évaluer la population totale d’élèves impactés.

Quelles informations exploiter ?

Le champ d’appellation officielle servira d’étiquette à l’établissement sur la carte.

Nous pourrons filtrer les établissements par leur état pour ne sélectionner que les établissements ouverts et utiliser le code nature de l’établissement pour mettre en évidence les écoles et collèges qui sont un ensemble de la population plus sensible à une exposition aux pesticides.

Finalement, notez que les coordonnées X et Y projetées sont fournies par l’IGN sur les données de décembre 2013. (Vous pouvez déjà remarquer dans le descriptif des données que plusieurs projections sont utilisées)

Données des parcelles de cultures en France métropolitaine

Où se procurer les données ?

Les Etats Membres de la Communauté Européenne ont l’obligation de localiser et d’identifier les parcelles agricoles, nous profiterons de l’opportunité pour télécharger le registre parcellaire 2012 mis à disposition par l’Agence de services et de paiement sur la plateforme data.gouv.fr.

Téléchargez toutes les ressources. Les fichiers sont distribués sous formes d’archives compressées à raison d’un fichier par département, ce qui vous demandera un peu de rigueur et beaucoup de patience !

Vérifiez que vous avez le bon nombre de départements, vous devriez avoir 96 ressources si vous téléchargez également la description des fichiers et les codes groupes cultures : la corse est scindée en départements 2A et 2B (il n’y a pas de numéro 20), le dernier numéro de département en France métropolitaine est le 95, le département 28 soumis en février 2014 est vide cependant les données sont présente dans la ressource intitulée « No Name », et les départements 75 et 92 ne possèdent pas de fichiers, même vides, ce qui doit refléter une absence de cultures déclarées. (Les quelques 1200 m2 de vignes parisiennes ne devraient pas avoir d’impact significatif sur les résultats ?!)

Les départements d’outre mer ne sont pas fournis, en l’absence de données je ne traiterai pas ces cas.

Quelles informations exploiter ?

Chaque fichier archive contient des données spécifiques aux îlots de culture pour le département indiqué par son code présent dans le nom d’archive.

Les informations de géométrie des îlots sont reliées à deux attributs:

  • Le numéro d’attribut Num_Ilot sert d’identifiant, il est unique
  • Le code groupe culture Code_Groupe_Culture est un attribut qui permet de déterminer le type de culture par rapport à la table de référence des codes groupes de culture. Cette table est fournie sous forme de fichier csv, elle permet d’associer un descriptif à la culture effectuée sur l’îlot et un code couleur. Le type de culture nous permettrait éventuellement d’associer aux données géographique un indicateur de fréquence des traitement (IFT) qui servirait à estimer les traitements phytosanitaires – mais pas la toxicité – pour chaque îlot. Les IFT seront récupérés de l’enquête Agreste Pratiques culturales 2011 et depuis une publication sur le site du Ministère de l’agriculture où figurent des données Agreste pour la viticulture en 2006 absentes sur l’enquête 2011

Données OpenStreetMap (OSM)

OpenStreetMap est un projet communautaire de grande envergure qui vise à créer et mettre librement à disposition et sous licence libre une carte du monde. OSM est alimenté par des milliers de contributeurs bénévoles qui maintiennent à jour les données cartographiques des réseaux routiers, ferroviaires, pédestres, les adresses, bâtiments et points d’intérêt sur la planète entière. OSM est exact et très à jour, il fournit des données que l’on peut exploiter pour nos propres besoins cartographiques.

Nous utiliserons les services d’OSM pour le rendu des fonds de carte dans l’application, et également au travers du logiciel QGIS pour la préparation du projet.

Données de l’Institut National de l’Information Géographique et Forestière (IGN)

L’IGN fournit sous licence ouverte la description de l’ensemble des unités administratives de France métropolitaine. Les données existent en 2 versions: GEOFLA® Communes, et GEOFLA® Départements, la réutilisation est gratuite pour tous les usages, nous allons exploiter ces données dans nos requêtes pour représenter l’information par département ou par commune.

Téléchargez GEOFLA® au format Shapefile, j’ai pris pour la réalisation les données 2013 des communes France métropolitaine et l’édition 2013 des départements France métropolitaine

Le logiciel QGIS Desktop

QGIS est un système d’information géographique (SIG) libre et open source qui permet de créer, éditer, visualiser, analyser ou encore publier des informations géographiques. Le logiciel QGIS Desktop est l’application pour ordinateur de bureau, disponible sur Windows, Mac, ou encore Linux. Cet outil, bien que gratuit, possède les capacités équivalentes à celles des logiciels propriétaires. Il est supporté par la fondation OSGeo et très largement utilisé.

Pour la réalisation de l’application cartographique je vais utiliser la version 2.4 du logiciel afin de préparer les couches géographiques qui seront exploitées sur le projet.

Je vous recommande l’excellent tutoriel QGIS créé par l’unité ADESS (Aménagement, Développement, Environnement, Santé et Société) sous tutelle du CNRS pour permettre aux débutants de d’initier aux SIG au travers de QGIS. Vous trouverez notamment en français des explications accessibles et illustrées sur la géodésie qui vous permettrons d’appréhender les systèmes de coordonnées et projections et qui pourrons vous servir de référence.

Utiliser un fond de carte OpenStreetMap

Nous allons voir comment l’installation d’une extension peut nous permettre de faire appel aux serveurs OSM afin de bénéficier d’un fond cartographique pour nos données.

Ouvrez QGIS, parcourez la barre de menu

Extension > Installer / Gérer les extensionsUne fenêtre s’ouvre, dans la zone de recherche saisissez le nom de l’extension OpenLayers Plugin

Installez l’extension, puis parcourez la barre de menuInternet > OpenLayers Plugin > OpenStreetMap > OpenStreetMap

Vous devriez désormais retrouver une couche nommée OpenStreetMap dans l’explorateur et la carte du monde doit progressivement apparaître dans la zone de rendu.

Le survol de la souris sur le code EPSG en bas à droite nous informe sur le système de coordonnées de référence, WGS 84 / Pseudo Mercator, et indique que le système de projection à la volée est activé.
Vous pouvez activer le système de projection à la volée dans les propriétés du projet (soit en passant par l’icône en bas à droite où apparaissent des ellipses, soit en passant par le menu Projet > Propriétés du Projet), la propriété s’active sous l’onglet SCR (système de coordonnées de références) et va permettre d’importer des couches sous un autre système de coordonnées.

Ces couches seront superposées dans le système de coordonnées défini par défaut pour le projet, ce qui nous sera utile pour l’import des couches dans le système de référence officiel pour la France métropolitaine. Les coordonnées des données d’établissements et d’îlots de culture sont fournies pour le système RGF93 – Lambert 93.

Ouverture des données d’établissements scolaires

QGIS propose d’ajouter des couches de texte délimité, ce qui permet de charger des fichiers au format texte ou texte csv à partir du moment où le fichier respecte un format de délimitation des champs et possède des champs de coordonnées géographiques.

Vous pouvez passer par le menu Couche > Ajouter une couche de texte délimité, ou bien utiliser l’icône en forme de virgule sur la palette à gauche, ou plus simplement faire glisser le fichier sur la fenêtre applicative de QGIS. Pour cette dernière option le logiciel essaie de déterminer lui même les paramètres, elle ne fonctionne pas dans notre cas.

Ajoutez une couche de texte délimité en passant par le menu, une fenêtre s’ouvre et vous permet de saisir des paramètres :

  • Choisissez le codage ISO-8859-1, il s’agit du jeu d’encodage des caractères Latin dans lequel se présente le fichier texte.
  • Le caractère délimiteur est le pipe ‘|’, vous devez décocher tous les délimiteurs prédéfinis et saisir le caractères dans la zone de saisies Autres délimiteurs(combinaison de touches Alt Gr + 6)
  • Précisez que la première ligne doit être ignorée en cochant l’option en-têtes en 1ère ligne
  • La définition de la géométrie est le Point, vous devez sélectionner les champs du fichier pour les valeurs de X et Y
Le logiciel vous averti de l’absence de géométrie sur 120 enregistrements… continuez.Au moment de sélectionner le système de coordonnées de référence, choisissez RGF93 / Lambert-93

Double cliquez sur la couche qui vient d’apparaître dans l’explorateur à gauche pour accéder aux propriétés. Sur la fenêtre qui s’ouvre, dans l’onglet Étiquettes choisissez d’étiqueter la couche avec le champ code_postal_uai, le rendu de carte va faire apparaître les codes des départements annexés aux points qui représentent les départements.

Abaissez le niveau de zoom de la carte afin de faire apparaître les établissements des départements d’outre mer (ç-a-d les codes département qui débutent par 97 ou 98). Que se passe-t-il ?

Les points qui correspondent aux départements d’outre mer n’apparaissent pas au bon endroit sur la carte, simplement parce que le fichier de données mixe les projections comme indiqué dans la documentation en annexe au fichier. Cela implique de traiter en amont le fichier.

Le chargement du fichier sous QGIS nous a fait prendre connaissance des problématiques des données d’établissement qu’il faudra prendre en compte lors du chargement du fichier en base de données. Premièrement certaines lignes ne possèdent pas de coordonnées géographiques, ce qui se vérifie aisément sous un tableur où l’on constate bien 120 enregistrements sans coordonnées. Deuxièmement il faudrait créer un import pour chaque système de référence employé, et malheureusement le seul champ qui nous permet cette discrimination est le code postal. Ultérieurement nous procéderons au chargement des données en base en restreignant les données à la France métropolitaine puisque nous ne possédons de toute façon pas les informations parcellaires des départements d’outre mer.

Ouverture des données de parcelles agricoles

Prenons l’archive RPG_2012_003.zip, le numéro nous indique qu’il s’agit de l’Allier. Si l’on décompresse l’archive, on dispose de 5 fichiers, dont 4 fichiers appartiennent au format SIG devenu standard : ESRI ShapeFile.

  • Le fichier Excel d’extension .xls nous renseigne sur les codes groupes culture comme nous avons pu voir lors de la récupération des données.
  • Le fichier d’extension .prj contient les informations sur le système de coordonnées, vous pouvez l’éditer avec un éditeur de texte pour obtenir des renseignements sur ce système
  • Le fichier .dbf stocke les données des attributs, vous pouvez l’ouvrir avec un tableur
  • Le fichier .shx stocke les indexes des enregistrements du fichier d’entités géographiques
  • Le fichiers .shp contient la géométrie des entités géographiques, il peut s’agir de points, lignes ou polygones. Dans ce cas précis il s’agit de contours de parcelles donc de polygones
Ouvrez QGIS ou si c’est déjà fait créez un nouveau projet.Faites glisser le fichier .shp sur la fenêtre d’application QGIS, le logiciel réalise l’opération d’ajout d’une couche vectorielle au projet.

Le résultat n’est pas en soit très parlant, nous pouvons ajouter la couche vectorielle des limites administratives départementales voir des limites communales afin d’obtenir une meilleure représentation.

Faites glisser les fichiers DEPARTEMENT.SHP, LIMITE_DEPARTEMENT.SHP, COMMUNE.SHP et LIMITE_COMMUNE.SHP que vous avez récupéré précédemment des données GEOFLA® sur la fenêtre d’application QGIS, vous pouvez contrôler les couches ajoutées dans l’explorateur à gauche et les réordonner pour faire apparaître les parcelles sur les autres couches.Vous pouvez contrôler les attributs en modifiant le style par type de culture.

Ouvrez les propriétés de la couche RPG_2012_003 en double cliquant sur le nom de la couche dans l’explorateur. Dans l’onglet Style de la fenêtre qui vient de s’ouvrir sélectionnez l’entrée Catégorisé dans la liste déroulante puis la colonne CULT_MAJ. Sélectionnez Couleurs au hasard pour la palette de couleurs puis cliquez sur le bouton Classer.

Vous devriez obtenir un résultat similaire à celui ci :

L’import des fichiers de données sous QGIS nous a permis de contrôler les formats, vérifier les attributs et détecter les anomalies. Nous n’allons pas directement exploiter les fichiers de données, passons maintenant au chargement des données en base de données spatiale pour bénéficier des possibilités d’indexation, d’accès concurrentiel et de requêtes spatiales.

Une fois les données en base nous retournerons sous QGIS avec pour objectif cette fois ci de créer un projet qui pourra être utilisé par un serveur cartographique pour transmettre des flux de données sur le réseau à destination de notre application de cartographie web.

Carte d’estimation d’exposition aux pesticides par établissement scolaire

Etude de la mise en oeuvre d’une application cartographique

Introduction

Ce projet a pour but d’étudier la mise en oeuvre d’une application cartographique.

J’ai choisi un sujet d’étude concret : il s’agira d’établir une carte d’estimation d’exposition aux pesticides pour les établissements scolaires en fonction de leur proximité aux zones agricoles.

Mes objectifs sont de documenter l’installation de l’architecture du système d’information et de construire une application de représentation de données qui démontre des cas d’utilisation statistiques courants.

Sommaire des points abordés

  • Définition d’un sujet d’étude
  • Base de données spatiales PostgreSQL PostGIS
  • Serveur cartographique QGIS Server
  • Application SIG QGIS Desktop
  • Cartes interactives sous Leaflet et OpenLayers 3
  • Requêtes spatiales

Aperçu du résultat :

La carte n’est pas accessible depuis Internet faute de moyens pour un hébergement suffisamment dimensionné.

Définition d’un sujet d’étude

Contexte

Exposition des enfants des zones rurales aux perturbateurs endocriniens

L’idée de l’application m’est apparue quand j’ai pris connaissance de cas d’intoxication aiguë d’enfants à l’école.

L’enquête EXPPERT (EXposition aux Pesticides PERTurbateurs endocriniens) menée par l’association Générations Futures sur les perturbateurs endocriniens s’intéresse plus particulièrement dans son 3ème volet à l’exposition des enfants qui vivent ou vont à l’école dans des zones agricoles.

Dans le cadre de l’enquête menée fin 2013, les prélèvements de mèches de cheveux ont été effectués sur 30 enfants et confiés pour analyse à un laboratoire de recherche indépendant.

La présence de perturbateurs endocriniens en nombre important dans les cheveux analysés révèle la forte exposition des enfants testés dans les 3 mois qui précèdent le prélèvement.

Les résultats peuvent être consultés dans le rapport d’étude.

Le rapport propose cette définition des perturbateurs endocriniens :

Les perturbateurs endocriniens (PE) sont des substances chimiques qui peuvent interférer avec le fonctionnement du système endocrinien et induire des effets néfastes sur l’organisme d’un individu ou sur ses descendants.

Les cancers hormonaux-dépendants (prostate, testicule, sein), les perturbations du métabolisme (obésité, diabète), de la reproduction (diminution de la fertilité, puberté précoce chez les filles), les problèmes cardiovasculaires mais aussi les troubles mentaux et du comportement, sont tous des effets potentiels des PE

Le rapport cible les enfants car ils font partie au même titre que les femmes enceintes des sous groupes de la population les plus sensibles, de la petite enfance jusqu’à la puberté.

Je fais le choix de ce sujet car il est d’actualité et préoccupant, parce qu’une représentation cartographique peu contribuer à faire prendre conscience du problème et apporter des informations à chacun pour mener ses propres réflexions.

Pour plus d’informations sur les effets des pesticides sur la santé vous pouvez consulter l’expertise de l’INSERM publiée en juin 2013.

Le projet de loi d’avenir pour l’Agriculture, l’Alimentation et la Forêt (LAAF)

Plusieurs associations ont déposé des propositions d’amendements pour la discussion de la LAAF en 2014.

Parmi les propositions on retrouve la volonté de protéger les populations riveraines des zones d’agriculture intensive par des zones de non traitement.

C’est le cas par exemple des propositions de l’AMLP (Alerte des médecins sur les pesticides).

La ministre de l’Ecologie en poste a proposé d’interdire toute pulvérisation de pesticides dans un rayon de deux cents mètres autour des écoles, suite à quoi la FNSEA (Fédération nationale des syndicats d’exploitants agricoles) aurait avancé des chiffres sur les millions d’hectares impactés par l’interdiction de l’usage des pesticides dans un rayon de 200 mètres de toutes les zones habitées.

Désormais, la LAAF conditionne l’autorisation de pulvérisation de pesticides près des écoles à des mesures de précaution : plantation de haies, gestion des horaires…

Ces actualités me font m’interroger sur les méthodes pour produire ces chiffres ou déterminer ces mesures, et je me pose la question de la prise en charge du problème par les ministres de l’agriculture et de l’écologie et l’absence de communication de la santé publique.

L’un des objectifs de l’application sera de permettre de se faire une idée personnelle sur la possibilité de répondre à ces interrogations sur les rayons de zones de non traitement et hectares potentiellement impactés en fonction de la distance de précaution.

Objectifs

A la lecture d’une étude de risques sanitaires publiée sur le portail InVS concernant les épandages de phytosanitaires menée en 2006, j’ai pu constater que la démarche adoptée consistait à recouper géographiquement les plaintes de la population avec les données de déclarations d’épandage sur parcelles.

Ce recoupement a été suivi d’une mesure de la qualité de l’air avec recherche de substances actives qui permet d’estimer la concentration dans l’air de ces substances hors période d’épandage et en période d’épandage. (A noter que les traitements aériens et traitements au sol étaient pris en compte en ce qui concerne cette étude).

Les concentrations dans l’air permettent d’estimer le niveau d’exposition sur la base de scénarios d’exposition aiguë ou chronique pour chaque valeur toxicologique de référence (VTR).

Cette étude suit à priori la conduite d’une évaluation des risques sanitaires telle que préconisée par le National Research Council américain en 1983, on peut reprocher à ce type d’études leur manque de recul sur les autres faits vérifiables et en cela l’initiative de Génération Futures apporte une approche moins théorique et met en avant la problématique par les faits.

La démarche de l’étude EXPPERT est pragmatique, elle s’appuie non pas sur des scénarios d’exposition mais sur des résultats d’analyse directement sur la population potentiellement exposée.

Bien que les résultats mettent à la lumière la problématique, l’étude est limitée par la sélection des échantillons sur la base du volontariat et le faible nombre analysés.

  • La carte que je propose d’établir prendra en compte uniquement l’exposition environnementale, il s’agira de mettre en évidence la proximité des établissements scolaires aux zones agricoles.
  • Comme l’age est un facteur de risque aggravant, l’application permettra de pointer plus particulièrement la petite école et les collèges.
  • Le type de culture influe également fortement sur les herbicides et insecticides déversés aussi la carte fera apparaître le type de culture, certaines cultures se prêtant d’avantage à des traitements toxiques, sans toutefois considérer les rotations de culture

Ces informations cartographiques pourrait servir à un échantillonnage statistique dans le cadre d’une étude d’impact plus significative et en mesure de représenter l’exposition des enfants français.

S’agissant des voies d’exposition aux pesticides les substances pénètrent selon soit la voie cutanée, soit la voie digestive soit la voie respiratoire : nous sommes tous exposés à un certain degré, mais pour la population générale la voie orale est considérée comme la voie d’exposition la plus importante.

L’enquête porte sur la proximité des enfants aux zones agricoles pulvérisées par des pesticides que ce soit à leur résidence principale ou à l’école.

J’aurai trouvé intéressant de constituer un échantillon d’enfants dont la résidence est éloignée des zones agricoles mais l’école à proximité.

Effectuer les mêmes analyses de cheveux sur les personnes au foyer des enfants concernés pourrait sans doute permettre de mettre en évidence l’exposition environnementale par voie respiratoire.

Cette démarche complémentaire permettrait d’avoir des informations sur la part d’exposition à la maison de celle de l’exposition en collectivité.

  • L’une des possibilités de l’application sera de donner la distance d’une adresse postale par rapport à une zone agricole.

Un autre objectif sera de déterminer par une requête spatiale le nombre d’hectares qui seraient impactés par une interdiction de l’usage des pesticides dans un rayon de 200 mètres des établissements scolaires.

Limitations du prototype

  • La carte ne couvrira que la France métropolitaine : les données de parcelles cultivées sont proposées sous une projection adaptée à la France métropolitaine, les départements et territoires d’outre mer doivent faire l’objet de données sous d’autres formats de projection dans des fichiers de données séparés
  • La précision de la carte est liée à la justesse des données. Pour les cultures il s’agit de déclarations, comme nous le verrons ultérieurement les données peuvent ne pas être à jour sur le type de culture ou encore absentes des cultures déclarées
  • Enfin il faut bien comprendre que la représentation propose des indications construites sur la base des informations disponibles de proximité géographique aux zones cultivées, il ne s’agit pas d’un indicateur de risque avéré par des mesures sur le terrain. Le manque de données sur les épandages par parcelles ne permet que des approximations et les risques réels d’exposition ne sont que des hypothèses