Tuilage des données OSM avec Leaflet

Carte d’occupation du sol (Landuse d’OSM)

osm landuse diouck

Principe du tuilage:

Le tuilage est le re-découpage d’un raster en petite images à des échelles prédéfinies pour être utilisé par des outils web.

Le tuilage va donc découper le raster sur plusieurs échelles (pyramidage) en autant de petites images (tuilage), dont le nom du dossier et celui de l’images sont les coordonnées x et y dans le système de projection mercator.

diouck pyramides

1/- Données: Pays de la Loire à télécharger depuis ce lien

2/- Analyse thématique uniquement sur la donnée « landuse »  avec Qgis et génération d’un fichier  au format png que vous pouvez télécharger  là :  OSM Landuse.

3/- Création de mosaïques d’images tuilées avec gdal2tiles.

Comme les autres outils de GDAL , il est assez facile à utiliser. Ci-dessous est un exemple de commande pour créer des tuiles de notre image.

gdal2tiles.py -z -v 7-13 "OSM Landuse.png" "OSM Landuse"

Pour faire simple nous allons créer une pyramide avec les niveaux de zoom allant de 7 à 13 dans le dossier OSM Landuse. La derniere version de gdal vous crée automatiquement un dossier pour chaque niveau de zoom.

4/- Intégration des tuile dans une application.

Je ne vais pas rentrer dans les détails de créations de cartes mais juste juste vous monter comment intégrer notre pyramide tuilée dans une application. Vous avez simplement à rajouter cette ligne:

var mytile =L.tileLayer(‘landuse/{z}/{x}/{y}.png’,{tms: ‘true’}).addTo(map);

La carte est accessible depuis ce lien

Publicités

#Monalisa2PostGIS: Importation de nombreux raster (Traitement par lots)

monalisa

Pour ceux qui travaillent avec Postgis de manière plus ou moins fréquente ,  si vous avez des difficultés d’intégration de plusieurs raster en même temps, voici un  script batch  qui fera office d’ETL. Évidemment les applications sont multiples par exemple regrouper plusieurs dalles d’ortho ou d’images satellites et calculer des indices de végétation.

Le scripte peut être exécuté sans problème malgré un nombre important de raster ( Tester sur un dossier de 50 dalles de 1km).  Voici les paramètres à modifier:

  1- D’abord définir les paramètres de connections

REM Parametre de connections
set PGPORT=5432 
set PGHOST=localhost
set PGUSER=postgres
set PGPASSWORD=postgres

2- Définir les  exécutables de postgres

REM  chercher le programme 
cd "C:\Program Files\PostgreSQL\9.5\bin"

3- Transformer en SQL l’image dans le schéma raster

Monalisa est disponible sur sa page  wikipedia

raster2pgsql   -d -I -C -e -Y   -s 4326   -t 128x128    "C:\raster\monalisa\*.jpg"  raster.monalisa >    "C:\raster\monalisa.sql" 
Avec -d pour supprimer la table si elle existe
-C pour la créer si elle n'existe pas
-I pour créer un overview
-s pour la projection et -t pour le tuilage

4- Importer le SQL dans la base de donnee bdu

psql -d bdu -f "C:\raster\monalisa.sql"

Pour info le script bat est disponible ici

5- Et visualiser votre résultat sur Qgis.

 

 

 

 

 

Calcul de l’indice de végétation NDVI avec PostgeSQL/Postgis

En télédétection, les indices font parties des méthodes de traitement que l’on appelle les transformations multispectrales. Ils consistent à convertir les luminances mesurées au niveau du capteur satellitaire en grandeurs ayant une signification dans le domaine de l’environnement.

Basés sur le caractère multispectral des données satellitaires, ils permettent de décrire l’état d’un phénomène. Un indice de végétation par exemple, peut rendre compte du stade de croissance végétale à un moment donné.
Tous les indices, que ce soient les indices de végétation, les indices des sols, les indices relatifs à la colonne d’eau, etc., reposent sur une approche empirique basée sur des données expérimentales. Les indices de végétation sont très utilisés d’une part, pour identifier et suivre la dynamique de la végétation, mais aussi pour estimer certains paramètres biophysiques caractéristiques des couverts végétaux, comme la biomasse, l’indice de surface foliaire, la fraction de rayonnement photosynthétique actif, etc.

L’indice le plus connu et le plus utilisé est l’indice de végétation par différence normalisé ou indice de Tucker (NDVI en anglais) (Rouse and Haas, 1973 ; Tucker, 1979). Son expression est la suivante :

(pir-rouge)/(pir+rouge)

Je suppose que vous vous êtes déjà familiarisés avec Postgis ,Qgis et Gdal.

Les données : Images Quickbird (Visible +proche infrarouge)

composition colorée DIOUCK

Composition colorée Rouge +PIR

Assemblage des données avec Gdal

Pour faciliter les traitements, j’ai assemblé les bandes du rouge et du proche infra-rouge dont nous aurons besoin pour calculer le NDVI. Si vous n’avez pasGdal je vous conseille d’intaller FWtools

–Stack Bande
gdal_merge -separate  bande3.tif bande4.tif -o bande34.tif

Pour vérifier notre stack voici la commande gdalinfo bande34.tif

gdalinfo_bande

 gdalinfo sur le stack

Et là on distingue bien nos deux bandes : La bande 1 correspond au rouge et la bande 2 au PIR

Importation des données dans Postgis :

J’ai créé un script bat pour l’intégration des données

REM Paramètres de connections
set PGPORT=5432
set PGHOST=localhost
set PGUSER=postgres
set PGPASSWORD=postgres
REM on traite le raster en le transformant en une suite de commandes SQL pour le schema ‘raster’ . Le –t permet de tuiler la donnée et le –I une indexation pour  gagner en performance
« C:\Program Files\PostgreSQL\9.3\bin\raster2pgsql »  -s 4326 -C -I -r -M -F -t 100×100  « D:\Postgis Raster\Satellite\QuickBird\ bande34.tif »  raster. bande34> bande34.sql
REM on crée la table raster dans notre base BDU de donnée via psql
« C:\Program Files\PostgreSQL\9.3\bin\psql » -d bdu -f bande34.sql
 
 
 
 

Calcul de l’indice de végétation  NDVI

L’indice le plus connu et le plus utilisé est l’indice de végétation par différence normalisé ou indice de Tucker (NDVI en anglais) (Rouse and Haas, 1973 ; Tucker, 1979). Son expression est la suivante :

La formule est la suivante (pir-rouge)/(pir+rouge)

Voici la requête dans postgis

create table ndvi as
select ST_MapAlgebra(a.rast, 2, a.rast, 1, ‘([rast1] – [rast2])/([rast1] + [rast2])::float’, ’32BF’) as rast
FROM raster.bande34 a

 

Visualisation dans QGIS

Comme vous pouvez le constater nous avons bien les valeurs du NDVI comprises entre -0.9 à 0.68.
Une simple classification permet de visualiser les différentes classes d’occupation du sol qui sont pareilles que la classification que j’ai fait avec les mêmes données à partir de QGIS

NDVI avec POSTGIS

NDVI POSTGIS DIOUCK
NDVI avec QGIS
NDVI QGIS DIOUCK

 

 

Pyramide avec GDAL et GEOSERVER

Retour d’expériences

Pyramide tuilé avec GDAL
Suite à un travail faite sur l’optimisation d’ortho photos de Rennes Métropole, je partage avec vous mon expérience sur la valorisation et la production d’images tuilées avec des overview et la mise en place d’une structure pyramidale d’image tuilées

Pour plus de documentation complète je vous renvoie à Geoserver Stéroide
Pour commencer télécharger gdal

Il y a un bug dans la version du script gdal_retile.py dans FWTOOLS 2.4.7.
Si vous rencontrer cette erreur
File "C:\PROGRA~1\FWTOOL~1.7\bin\gdal_retile.py", line 273
print("Building internal Index for %d tile(s) ..." % len(inputTiles), end=' ')
^
SyntaxError: invalid syntax

Remplacer la ligne 273 par
print "Building internal Index for %d tile(s) ..." % len(inputTiles)
Test GDAL avec FWTOOLS : 1329_7218.tif
Je ferais mes tests à partir une dalle geotiff : 1329_7218.tif
Dans un premier temps nous allons consulter l’information de base contenue dans notre image.
Pour ce faire nous allons utiliser la commande :
Gdalinfo 1329_7218.tif

gdal_info

Comme vous pouvez le constater nous avons en gros une Geotiff avec une projection en Lambert93, des block de 5000*5000

Pour optimiser les données, la bibliographie conseil d’utiliser des des blocks de 256*256 avec des overview internes ou externes.

Pour créer des blocks vous devez utiliser Gdal_translate qui est l’outil de transformation de GDAL (to translate = transformer, convertir), il sert principalement à changer le format, mais aussi à extraire une dalle à partir d’un raster existant.
gdal_translate -co "TILED=YES" -co "BLOCKXSIZE=256" -co "BLOCKYSIZE=256" "C:\Documents and Settings\a.diouck\Bureau\GDAL\ortho_15cm\assemblage\beton\depart\20cm.tif" "C:\Documents and Settings\a.diouck\Bureau\GDAL\ortho_15cm\assemblage\beton\depart\20cmblock256.tif"
gdaltranslate

Pour vérifier la transformation vous faites

Gdalinfo 1329_7218.tif

gdalinfo2

Maintenant avec notre image nous avons à la fois des block de 256*256 mais aussi des overview internes que l’on peut directement utiliser dans un serveur cartographique pour créer des tuiles performantes.
Mais avant cela je vais vous montrer quelque chose de pratique c’est-à-dire la reprojection des données avec la commande gdalwarp, utile si on veut changer de coordonnées.
Gdalwarp est l’outil de « déformation » de GDAL (to warp = déformer), il sert principalement à reprojeter, mais aussi à modifier les pixels de l’image.

* -of pour préciser le format (ex : -of ECW)
* -co pour appeler une option de création et en préciser la valeur (ex : -co TARGET=90)
* -projwin dans gdal_translate, pour extraire une dalle d’un raster (ex : -projwin 575000 325000 576000 322000)

* -s_srs dans gdalwarp, pour spécifier la projection du fichier source (ex: -s_srs EPSG:27562)
* -t_srs dans gdalwarp, pour spécifier la projection du fichier en sortie (ex: -s_srs EPSG:3948)

Example passage de rgf 93 par defaut de notre image source vers le CC48 utilisée par Rennes Métropole
gdalwarp -t_srs EPSG:3948 1329_7218.tif 1329_7218_3948.tif
gdalinfo 3

Visualisation avec Qgis
qgis gdal

Structure pyramidal
Maintenant que vous savez comment manipuler gdal pour optimizer vos tuiles, je vais vous montrer une autre facon plus traditionnelle de creation de tuiles .
Gdal nous permet de créer des pyramides d’apres mon expériance encore plus performant que ce qu’on vient de voir.
Pour ce faire j’ai regroupé les commande de tout à l’heure plus gdal_retile pour générer une pyramide avec une reéchantillonnage bilineaire, sur 4 niveau, des dalles de 2048*2048 , des block de 256*256 dans un dossier bmpyramid
gdal_retile.py -v -r bilinear -levels 4 -ps 2048 2048 -co « TILED=YES » -co « COMPRESS=JPEG » -co « BLOCKXSIZE=256 » -co « BLOCKYSIZE=256 » -targetDir bmpyramid 1329_7218_3948.tif

gdal_retile

Comme vous pouvez le constater nous avons 4 niveaux et le reste c est le niveau “0”
pyramide

Construction d’un pyramide avec Geoserver
Un autre example

gdal_retile -v -ps 500 500 -s_srs EPSG:3948 -levels 7 -tileIndex 20cm_optimise.shp -tileIndexField tilename -targetDir « C:\FWTools2.4.7\geotif\pyramide\arrivee » « C:\FWTools2.4.7\geotif\pyramide\depart\20cm.tif »
gdal_retile 2
Tous les images en dehors des dossiers sont le niveau « 0 » . En important le dossier depuis Geoserver avec « imagepyramid » Geoserver va créer un dossier “0” pour regrouper ltous les images en dehors des dossiers.pyramide

pyramide 2

Il crée également des indexes spatiales et fichiers de projection pour chaque dossier.
Et enfin à la racine il crée un fichier de propriété qui fera le lien entre chaque dossier.
geoserver