Archives par mot-clé : Qgis

Un conteneur pour postgreSQL et postgis

Il arrive des moments dans la vie où on a besoin d’un petit serveur de base de données rapidement, là tout de suite. Si ça vous est déjà arrivé, vous savez que c’est parfois se lancer dans une aventure … qui si elle n’est pas périlleuse, peut s’avérer longue et douloureuse.

Alors, pourquoi ne pas en profiter pour passer à docker ? Parce qu’avec le principe de dockerfiles, l’installation ne prend pas plus longtemps que de télécharger les packages.

Installation docker sur Ubuntu

Source : docker.com

Le site de docker est très bien fait. Sur Ubuntu 16.04, il suffira de lancer les lignes de commandes suivantes.

sudo apt-get install -y --no-install-recommends apt-transport-https ca-certificates curl software-properties-common
 curl -fsSL https://apt.dockerproject.org/gpg | sudo apt-key add - apt-key fingerprint 58118E89F3A912897C070ADBF76221572C52609D
 sudo add-apt-repository "deb https://apt.dockerproject.org/repo/ \
 ubuntu-$(lsb_release -cs) \
 main"
 sudo apt-get update
 sudo apt-get -y install docker-ce

Si tout s’est déroulé sans accrocs, vous pouvez vérifier que docker fonctionne avec :

sudo docker run hello-world

Il ne reste enfin quelques petits ajustements à faire, notamment configurer votre utilisateur dans le groupe docker pour s’abstraire de l’utilisation de sudo par exemple. Vous trouverez ces choses-là sur cette page.

Installation de postGIS et postgreSQL

Source : postgis dockerfile

Pour créer un conteneur docker il suffit de lancer la commande docker run. Docker va vérifier si vous disposez d’une « image » du service que vous souhaitez lancer. Si c’est le cas, un conteneur sera créé, sinon les différentes composantes de l’image seront téléchargées et le conteneur sera créé dans la foulée.

On commence ici par créer un conteneur postGIS qui s’appellera psql-futurSahel (C’est l’étiquette du conteneur et cela n’aura pas d’influence sur la base de données). On définit le mot de passe de postgreSQL et on spécifie que l’on veut utiliser l’image postgis fournie par l’utilisateur mdillon.

docker run --name psql-futurSahel -e POSTGRES_PASSWORD=postgres -d -p 5432:5432 mdillon/postgis

Pour que cela fonction, il faut également installer un postgreSQL et le lier à l’image de postGIS. Ce qui nous permettra ensuite d’accéder à l’interface de base de données et de nous connecter avec psql

docker run -it --link psql-futurSahel:postgres --rm postgres \
 sh -c 'exec psql -h "$POSTGRES_PORT_5432_TCP_ADDR" -p "$POSTGRES_PORT_5432_TCP_PORT" -U postgres'

Dans Postgre

Si tout s’est passé comme il faut, suite à la dernière commande, et après téléchargement des différents constituants de l’image, vous êtes rentré dans le conteneur et vous avez en face de vous un prompt psql.

Vous pouvez donc utiliser les commandes habituelles de postgre. Il faut, une fois connecté à psql, créer une base de données

CREATE DATABASE fs_gis;

On peut ensuite se connecter à la base de données

\connect fs_gis;

Et lancer la commande

CREATE EXTENSION postgis;

La base est donc créée en utilisant le template postGIS.

Peupler la base de données

Vous avez maintenant un serveur de base de données spatiales qui tourne et auquel vous pouvez accéder depuis la machine hôte.

Pour connaître l’IP du conteneur :

docker inspect idconteneur

Vous pouvez alors configurer vos outils de gestion de base de données en utilisant l’IP du conteneur et/ou votre propre IP (ifconfig).

Dans Qgis vous pouvez configurer la gestion de la base de la manière suivante :

Et vous pouvez utiliser l’extension Base de données > Gestionnaire de base de données > Gestionnaire de base de données, pour importer dans postgreSQL vos données en cliquant sur ‘importer une couche ou un fichier’. Cet outil va utiliser le script shp2psql de manière transparente pour l’utilisateur.

Earthdata des données SIG pour l’environnement

Voilà plusieurs mois que je n’avais pas pris le temps de faire un post sur mon blog !
On trouve toujours beaucoup d’excuses, des papiers sur le feu, un déménagement, un nouveau post-doctorat… Alors pour recommencer doucement, une petite note sur un service de la NASA que j’ai découvert aujourd’hui : Reverb | echo[1]

logo earthdata
Un peu de contexte à cette découverte, mon nouveau post-doc m’a amené à sortir de ma zone de confort, on s’éloigne de la viticulture et de l’eau (encore que), pour travailler sur la zone sahélienne du Sénégal. Je cherchais des indicateurs de désertification… Et je suis tombé sur un blog génial de Martin Brandt[2] (chercheur au Department of Geosciences and Natural Resource Management à l’université de Copenhagen (Denmark)). En plus de pointer sur des papiers scientifiques, il y a aussi quelques posts sur des outils de télédétection (GRASS-GIS) et des liens vers des pourvoyeurs de données. C’est là que j’ai découvert LE truc.

Reverb est un silo de données environnementales propulsées par la NASA qui contient des choses incroyables et un puissant moteur de recherche. Je vous engage donc à regarder cette petite vidéo d’explication[3] (en anglais).

Je ne reviendrai pas sur la manière de trouver ses données, la vidéo est là, par contre une chose géniale qui pourrait inspirer d’autres sites : une fois que j’ai fait mon petit marché, le service me propose de télécharger un fichier TXT contenant les URLs des ressources demandées. L’idée est génial. Il y a quelques semaines, j’avais du faire une centaine de copier-coller depuis le site opendata.gouv.fr pour pouvoir automatiser un traitement sur le RGP de 2012[4] (j’en parlerai peut-être dans quelques mois quand l’article sera publié).

Donc une grande satisfaction, je disais, pour cette bonne pratique… je voulais donc me lancer dans un script bash et du wget… mais les ennuis arrivent …

wget http://e4ftl01.cr.usgs.gov//MODV6_Cmp_B/MOLT/MOD09Q1.006/2016.09.05/MOD09Q1.A2016249.h16v07.006.2016258071050.hdf

Connexion à urs.earthdata.nasa.gov (urs.earthdata.nasa.gov)|198.118.243.33|:443… connecté.
requête HTTP transmise, en attente de la réponse… 401 Unauthorized
Échec d’authentification par identifiant et mot de passe.
URL transformed to HTTPS due to an HSTS policy

Mince un login ! Nooon ! Pourtant Martin Brandt ne semble pas avoir ce problème. La solution est en deux étapes. Dans un premier temps, il faudra se créer un compte « earthdata » en cliquant sur login en haut à gauche. Ensuite, et la solution est sur le wiki [5], il faudra faire un petit fichier de config « curl ».

cd ~
touch .netrc
echo "machine urs.earthdata.nasa.gov login ; password " > .netrc
chmod 0600 .netrc

Vous l’aurez compris, il ne reste qu’à remplacer par votre login et par votre mot de passe. Et cette fois-ci les requêtes passeront !

wget http://e4ftl01.cr.usgs.gov//MODV6_Cmp_B/MOLT/MOD09Q1.006/2016.09.05/MOD09Q1.A2016249.h16v07.006.2016258071050.hdf

Sauvegarde en : « MOD09Q1.A2016249.h16v07.006.2016258071050.hdf »
MOD09Q1.A2016249.h1 100%[===================>]  64,76M   344KB/s    in 1m 51s

Merci pour les données, et merci pour cette facilité d’utilisation !!!

Données MODIS issu de earthdata
Les liens

[1]https://reverb.echo.nasa.gov
[2]https://matinbrandt.wordpress.com
[3]https://www.youtube.com/watch?v=iBXzzUv3b4w
[4]https://www.data.gouv.fr/fr/datasets/registre-parcellaire-graphique-2012-contours-des-ilots-culturaux-et-leur-groupe-de-cultures-majorita/
[5]https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+cURL+And+Wget

pgRouting, postGIS, postgreSQL sont sur un bateau !

C’est un sujet que j’avais déjà traité dans un post antérieur, mais comme le blog de l’époque a disparu, il a fallu que je me replonge dans les bas-fonds de postgresql, postgis et pgRouting. Heureusement que les articles de Anita Graser (alias underdarck) sont encore en ligne, eux, mais ils ont nécessité quelques ajustements avec pgRouting 2!

Si vous en êtes là c’est que vous cherchez un moyen de faire des calculs d’isochrones et que, petit à petit, la solution pgRouting (même si elle peut sembler disproportionnée) s’impose à vous.

A ce stade, je partirais du principe que vous avez un serveur postgreSQL/postGIS fonctionnel. Vous trouverez plusieurs ressources sur Internet. Pour ma part, j’ai réalisé ce petit tuto sur Fedora 20 (un bon tuto pour l’installation postgresql postgis su le blog de Marcus Neteler et pgRouting est dans les packages de la distribution).

Attention toutefois à la version. Il se trouve que, sur Fedora, pgRouting est encore à sa version 1.05 qui nécessite quelques ajustement à la mano (révélé par Gregor the map guy). De mon côté, j’ai plutôt opté pour une compilation depuis les sources qui sont disponibles sur gitHub.

Créer un réseau « routable »

La première chose à faire est de préparer la base de données :

createdb --encoding=UTF8 --owner=delaye geodb
psql -d geodb -c "CREATE EXTENSION postgis;"
psql -d geodb -c "CREATE EXTENSION postgis_topology;"
psql -d geodb -c "CREATE EXTENSION pgrouting;"

Ensuite,  il faut importer un shapefile dans postgreSQL. Vous pouvez utiliser un client graphique comme shp2pgsql-gui ou deux petites lignes de commande (des infos ici )

#transformer en requête SQL un shapefile
shp2pgsql -I -t 2D -s 2154 reseau_total_AOC.shp acces_banyuls access_banyuls > acces_banyuls.sql
#envoyer les données dans la base 
psql -f acces_banyuls.sql -d geodb

Vous pouvez dès à present vous connecter à la base avec QGIS et avoir accès aux données

connection à postgreSQL avec QGIS 2.4
connection à postgreSQL avec QGIS 2.4

Pour pouvoir créer un graph connecté, il faut identifier pour chaque ligne ou polyligne le début et la fin. C’est ce que nous allons faire avec :

--evaluer  le type de la colonne geom
select * from geometry_columns;
--pour évaluer  la dimension des coordonnées de vos géométries
select distinct st_coordDim(geom) from acces_banyuls;

--convertir la colonne en deux dimensions si les données sont en 2D si vous n'avez pas utilisé l'option -t 2D de shp2pgsql
alter table acces_banyuls alter column geom type geometry(Multilinestring, 2154) using st_force2D(geom);

On va ensuite devoir transformer des multiLine-string en line-string dans postGIS parce que dans sa version 2 pgRouting ne travaille plus qu’avec ça .

--d'utilisation de st_dump pour transformer multiLine-string en line-string
CREATE TABLE simple_ways AS SELECT *, (ST_Dump(geom)).geom AS simple_geom FROM acces_banyuls WHERE geom IS NOT NULL;
--au passage on suprime la colonne multiline string
ALTER TABLE simple_ways DROP COLUMN geom;

On va ensuite créer une colonne temps de trajet à partir d’une colonne existante dans la base de données qui contient, pour chaque tronçon, la vitesse limite :

--mise à jour des temps de trajet
ALTER TABLE simple_ways ADD COLUMN length_m INTEGER;
UPDATE simple_ways SET length_m = st_length(st_transform(simple_geom,2154));

ALTER TABLE simple_ways ADD COLUMN traveltime_min DOUBLE PRECISION;
UPDATE simple_ways set traveltime_min = length_m / (vtss_pt / 60 * 1000)  --pour avoir la vitesse en metre

On peut ensuite essayer les fonctions de pgRouting proposer par Anita Graser

-- ajouter des col "source" et "target" à la table
ALTER TABLE acces_banyuls ADD COLUMN "source" integer;
ALTER TABLE acces_banyuls ADD COLUMN "target" integer;

--Creation de colonne startPoint et endpoint (syntaxe à changer ne pas oublier st_)
CREATE OR REPLACE VIEW road_ext AS
   SELECT *, st_startpoint(simple_geom), st_endpoint(simple_geom)
   FROM simple_ways;

--Une table qui contient tout les noeuds
CREATE TABLE node AS
   SELECT row_number() OVER (ORDER BY foo.p)::integer AS id,
          foo.p AS the_geom
   FROM (        
      SELECT DISTINCT road_ext.st_startpoint AS p FROM road_ext
      UNION
      SELECT DISTINCT road_ext.st_endpoint AS p FROM road_ext
   ) foo
   GROUP BY foo.p;

   
--creation du réseau  -> 128330 ms 
CREATE TABLE network AS
   SELECT a.*, b.id as start_id, c.id as end_id
   FROM road_ext AS a
      JOIN node AS b ON a.st_startpoint = b.the_geom
      JOIN node AS c ON a.st_endpoint = c.the_geom;

Pour la dernière commande il faut faire attention, en effet à chaque fois que vous regénérez une table node, les id de ces derniers peuvent changer. Donc quand vous allez devoir spécifier l’ID du noeud source (6791 dans la requête suivante), faites attention de bien le mettre à jour.

--Creation de la table qui contient les temps de trajet
--besoin d'une colonne avec travelTime pour faire le calcul
-- temps de calcul 123190 ms
create table cost_road as
select
    id,
    the_geom,
    (select sum(cost) from (
       SELECT * FROM pgr_dijkstra('
       SELECT gid AS id,
          start_id::int4 AS source,
          end_id::int4 AS target,
          traveltime_min::float8 AS cost
       FROM network',
       6791,
       id,
       false,
       false)) as foo ) as cost
from node;

Vous pouvez ensuite charger ces points dans Qgis et lancer une interpolation avec l’extension interpolation, par exemple, qui va vous en proposer deux types :

  1. TIN (Triangulated Irregular Netwaork)
  2. IDW (Inverse Distance Weighted)

Voilà le résultat d’une interpolation IDW (avec une discrétisation moche je vous l’accorde 🙂 mais on a qu’à dire que c’est pas l’objet ici).

isochrones

Ou, si vous voulez une interpolation plus propre, vous pouvez aussi passer sous R.