Table des matières
PAGE OBSOLETE

!!NB!! Depuis que GDAL est compilé avec PostGIS, l'outil raster2pgsql n'est plus un script python script. Le script python peut ne plus fonctionner avec les futurs version de PostGIS Raster (cf http://www.postgis.org/documentation/manual-2.0/using_raster.xml.html). Certaines fonctions ont par ailleurs été renommée dans la version finale.

, , ,

Installation de PostGIS-Raster/PostGIS 2.0 pour Windows XP

Installation de PostGIS 2.0

  1. Télécharger les binaires compatible avec la version de PostreSQL (http://www.postgis.org/download/windows/experimental.php)
  2. Lire les consignes des fichiers textes README et README_Raster:
    • Éditer makepostgisdb.bat
    • Préciser le mot de passe de connexion
    • Décommenter la ligne pour les templates, si nécessaire
  3. Lancer makepostgisdb.bat
  4. La base de données template_postgis20 est créée, avec le schéma public contenant les fonctions raster et la schéma topology

Installation de Python 2.7

http://www.python.org/

Installation du package NumPy et de la librairie GDAL (tentative 1)

http://www.lfd.uci.edu/~gohlke/pythonlibs Télécharger les versions compatibles et lancer les.

Tester l'installation en testant dans la fenêtre de commande raster2pgsql.py (cf NB tentative 2)

NB: suite à l'exécution de ces 2 packages, et ce moment là de l'installation, je n'exécutais vainement que raster2pgsql. J'ai donc testé une autre installation de GDAL. (cf ci-après)

Installation de la librairie GDAL (tentative 2)

via l'installateur générique de GDAL (http://vbkto.dyndns.org/sdk/Download.aspx?file=release-1310-gdal-1-8-mapserver-5-6\gdal-18-1310-core.msi)

Bien suivre l'étape 2 de la page http://gis4free.wordpress.com/category/postgis-raster/, notamment l'UPDATE 2011-03-25 sur l'ajout de variables d'environnement.

Tester l'installation en testant dans la fenêtre de commande raster2pgsql.py.

NB: au premier test, comme lors de la tentative 1, c'est raster2pgsql qui était vainement lancé, jusqu'à ce que raster2pgsql.py soit testé avec succès. A voir si lors de la tentative 1 raster2pgsql.py fonctionnait.

Notes sur certaines fonctions de PostGIS Raster (PostGIS 2.0)

http://www.postgis.org/documentation/

Accessoires Raster

ST_GeoReference, ST_Height, ST_MetaData, ST_ScaleX, ST_ScaleY, ST_SkewX, ST_SkewY, ST_SRID, ST_UpperLeftX, ST_UpperLeftY, ST_Width

(cf. http://en.wikipedia.org/wiki/World_file)

Les métadonnées de Raster utilisent 6 types de descripteurs :

ST_Raster2WorldCoordX, ST_Raster2WorldCoordY

Renvoie les coordonnées x/y, dans le référentiel spatial du raster, du coin haut gauche d'un raster, d'une colonne ou d'une ligne.

ST_World2RasterCoordX, ST_World2RasterCoordY

Renvoie la colonne d'un ponctuel (géométrie ou coordonnées) du même référentiel spatial que le raster.

ST_Histogram

Permets la classification des valeurs d'une bande donnée, et d'afficher de fait les valeurs minimum et maximum en excluant les valeurs '-1'. 

SELECT band, (stats).*
FROM (SELECT rid, band, ST_Histogram(ST_SetBandNoDataValue(raster,-1), band, 5) As stats
FROM bd_topo.mnt CROSS JOIN generate_series(1,1) As band
WHERE rid=1) As foo

-----resultats-----
'band'|'min'  |'max' |'count' |'percent'
1     |0      |253.6 |11004566|0.585324
1     |253.6  |507.2 |2861291 |0.15219
1     |507.2  |760.8 |3102829 |0.165037
1     |760.8  |1014.4|1667889 |0.088714
1     |1014.4 |1268  |164249  |0.008736

Temps d'exécution : 6 secondes

Accessoires et Manipulateurs de Pixels

ST_Value

Exemple ci-après pour l'énumération des valeurs des 100 premières lignes et colonnes , soit près de 5/10 000ème du MNT du Département de l'Hérault 

SELECT x, y, ST_Value(raster, 1, x, y) AS b1val
FROM bd_topo.mnt 
	CROSS JOIN generate_series(1,100) AS x 
	CROSS JOIN generate_series(1,100) AS y
WHERE rid = 1;

Temps d'exécution : 46 minutes

CAS D'UTILISATIONS

Caractéristiques des images utilisées:

SRTM (WGS84)

MNT de département de l'Hérault(Lambert93)

Vectorisation d'une sélection sur la valeur de pixels dans une sous partie du raster

Soit: vectorisation des altitudes de plus de 300 mètres dans une zone de 250 mètres de côté dans le nord-est de l'emprise du MNT du Département de l'Hérault (10 colonnes/lignes de 25 mètres depuis les coin haut-gauche), soit . 

 SELECT rid, ST_Union(pixpolyg)  As geometrie
FROM (
	SELECT rid,
		ST_Translate(
			ST_MakeEnvelope(
				ST_UpperLeftX(raster), ST_UpperLeftY(raster),
				ST_UpperLeftX(raster) + ST_ScaleX(raster),
				ST_UpperLeftY(raster) + ST_ScaleY(raster), 2154
				), 
			ST_ScaleX(raster)*x-ST_ScaleX(raster), 
			ST_ScaleY(raster)*y-ST_ScaleY(raster)
		) As pixpolyg, 
		ST_Value(raster, 1, x, y) As b1val
	FROM bd_topo.mnt 
		CROSS JOIN generate_series(1,10) As x 
                CROSS JOIN generate_series(1,10) As y
	WHERE rid = 1
	) As foo
WHERE b1val > 300::double precision
  GROUP BY foo.rid

Temps d'exécution : 4 minutes (moins de 3 minutes le lendemain)

NB: Dans la documentation, les exemples de ce type présentent la fonction ST_Translate avec les options de translation suivante 

ST_Translate(
			ST_MakeEnvelope(
				ST_UpperLeftX(raster), ST_UpperLeftY(raster),
				ST_UpperLeftX(raster) + ST_ScaleX(raster),
				ST_UpperLeftY(raster) + ST_ScaleY(raster), 2154
				), 
			ST_ScaleX(raster)*x, 
			ST_ScaleY(raster)*y
		)

causant un décalage d'un pixel en x et y!!!

création d'une couche vecteur issue de l'intersection d'un raster (récupération de la valeur du pixel) et d'un vecteur (récupération du gid)

SELECT foo.rid, foo.id_entite_occsol, ST_AsText((foo.geomval).geom) AS geomwkt, (foo.geomval).val
FROM(
	SELECT A.rid, ocsol_1_commune34.id_entite_occsol , ST_Intersection(A.raster, ocsol_1_commune34.geometrie) AS geomval
	FROM bd_topo.mnt AS A 
		CROSS JOIN dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text, '
			SELECT id_entite_occsol, area, perimeter, idpol_06, st_transform(raster_test_ocsol_1_commune34.geometrie, 2154) as geometrie, id_poste, code_insee
			FROM raster_test_ocsol_1_commune34', true
		) ocsol_1_commune34(id_entite_occsol integer, area double precision, perimeter double precision, 
			idpol_06 character varying, geometrie geometry, id_poste character varying, 
			code_insee character varying)
	WHERE A.rid =1 
) AS foo;

ou 

SELECT foo.rid, foo.id_entite_occsol, ST_AsText((foo.geomval).geom) AS geomwkt, (foo.geomval).val
FROM(
	SELECT A.rid, ocsol_1_commune34.id_entite_occsol , ST_Intersection(A.raster, ocsol_1_commune34.geometrie) AS geomval
	FROM bd_topo.mnt AS A 
		JOIN dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text, '
			SELECT id_entite_occsol, area, perimeter, idpol_06, st_transform(raster_test_ocsol_1_commune34.geometrie, 2154) as geometrie, id_poste, code_insee
			FROM raster_test_ocsol_1_commune34', true
		) ocsol_1_commune34(id_entite_occsol integer, area double precision, perimeter double precision, 
			idpol_06 character varying, geometrie geometry, id_poste character varying, 
			code_insee character varying)
			ON st_intersects(ocsol_1_commune34.geometrie, A.raster)
	WHERE A.rid =1 
) AS foo;

Ces deux requêtes font planter PostgreSQL et seul un redémarrage de l'ordinatuer permet la reconnexion au serveur!!!

Quant à 

SELECT id_entite_occsol, (gv).geom, (gv).val alti
FROM (
	SELECT id_entite_occsol, ST_Intersection(geometrie, raster) gv
	FROM dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text, '
			SELECT id_entite_occsol, area, perimeter, 
			idpol_06, 
			cleangeometry(st_transform(raster_test_ocsol_1_commune34.geometrie, 2154)) as geometrie, 
			id_poste, code_insee
			FROM raster_test_ocsol_1_commune34', true
		) ocsol_1_commune34(id_entite_occsol integer, area double precision, perimeter double precision, 
			idpol_06 character varying, 
			geometrie geometry, 
			id_poste character varying, code_insee character varying), 
		bd_topo.mnt
WHERE ST_Intersects(raster, geometrie)) sub
LIMIT 10

cela renvoie l'erreur 

ERREUR: Error creating GEOS Coordinate Sequence
État SQL :XX000
Contexte : PL/pgSQL function "_st_intersects" line 28 at affectation

ou un “Microsoft Visual C++ Runtime Library Error”

Test 2ème

Les mêmes opérations sont recommencées après que le raster ait été réintégré mais sous forme de tuiles de 250*250px.

Là les dernières requêtes aboutissent mais ne renvoient aucun résultat

Tests 3ème

SUCCES

Croisement de données de type ponctuel avec le MNT du département de l'Hérault : 

SELECT id_obs, nom_complet, cd_nom, 
		(gv).val, 
		(gv).geom AS geometrie
	FROM (
		SELECT id_obs, nom_complet, cd_nom, 
			ST_Intersection(raster, geometrie) AS gv
		FROM srtm.srtm_lr_wgs84_block,
			dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text,
				'SELECT id_obs, nom_complet, cd_nom, 
					st_transform(st_buffer(st_transform(geometrie, 2154), 100), 4326) as geometrie
				FROM saisie.saisie_observation
				LIMIT 10')
				tampon_pts_test(id_obs integer, nom_complet character varying, 
					cd_nom character varying, geometrie geometry)
		WHERE ST_Intersects(raster, geometrie)
		) foo

Résultat en 3.5 s: 

id_obs  |nom_complet                                            |cd_nom |val
2       |Zerynthia polyxena (Denis & Schiffermüller, 1775)      |8267	|518
3	|Pyrgus malvae (Linnaeus, 1758)	                        |53221	|59
28	|Tircis	Pararge aegeria (Linnaeus, 1758)	        |53595	|63
17235	|Issoria lathonia (Linnaeus, 1758)	                |53908	|1547
21248	|Chazara briseis (Linnaeus, 1764)	                |53425	|92
21690	|Gonepteryx cleopatra (Linnaeus, 1767)	                |54419	|233
21733	|Papilio sibilla Linnaeus, 1767	                        |392395	|391
22122	|Maculinea arion Linnaeus, 1758	                        |54085	|236
22123	|Maculinea alcon rebeli Hirschke, 1904	                |54083	|236
43	|Melanargia lachesis (Hübner, 1790)	                |219805	|142

NB: les données vecteur étant en WGS84, on réalise la zone tampon après avoir les avoir convertis dans un sytème projeté (ici Lambert93). PostGIS ne gérant pas la projection à la volée des raster, les zones tampons sont reprojetées en WGS84 pour réaliser le croisement spatial 

SELECT id_obs, nom_complet, cd_nom, 
		(gv).val, 
		(gv).geom AS geometrie
	FROM (
		SELECT id_obs, nom_complet, cd_nom, 
			ST_Intersection(raster, geometrie) AS gv
		FROM srtm.srtm_lr_wgs84_block,
			dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text,
				'SELECT id_obs, nom_complet, cd_nom, 
					st_transform(st_buffer(st_transform(geometrie, 2154), 100), 4326) as geometrie
				FROM saisie.saisie_observation
				LIMIT 10')
				tampon_pts_test(id_obs integer, nom_complet character varying, 
					cd_nom character varying, geometrie geometry)
		WHERE ST_Intersects(raster, geometrie)
		) foo
ORDER BY nom_complet

Temps d'éxécution : 3.5 s 

 premières lignes de résultats...
id_obs  |nom_complet                            |cd_nom |val
21248   |Chazara briseis (Linnaeus, 1764)	|53425	|70
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|71
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|78
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|79
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|92
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|102
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|109
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|126
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|148
21248	|Chazara briseis (Linnaeus, 1764)	|53425	|157
21690	|Gonepteryx cleopatra (Linnaeus, 1767)	|54419	|208
21690	|Gonepteryx cleopatra (Linnaeus, 1767)	|54419	|215
21690	|Gonepteryx cleopatra (Linnaeus, 1767)	|54419	|217
...     |...                                    |...    |...

Vectorisation du SRTM intersectant le département de l'hérault pour les altitudes >= à 600 m

SUCCES 

SELECT val, geometrie AS geometrie
FROM	(
	SELECT rid, (ST_DumpAsPolygons(raster)).geom AS geometrie, (ST_DumpAsPolygons(raster)).val
	FROM srtm.srtm_lr_wgs84_block
		JOIN dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text, 
				'SELECT code_dept, st_transform((geometrie), 4326) as geometrie
				FROM ign_bd_topo.departements_lr 
				WHERE code_dept=''34'''
				) dept_34(code_dept character varying, geometrie geometry)
		ON ST_Intersects(raster, geometrie)
	) srtm_34
WHERE val>=600

temps d'exécution : 5 minutes

un peu plus loin: Création d'une couche "forêts montagnardes" par croisement d'OCSOL et du SRTM
SELECT id_entite_occsol, id_poste, val, ST_Intersection(zone_montagne_34.geometrie, ocsol_foret_34.geometrie) AS geometrie
FROM	(
	SELECT val, geometrie
	FROM	(
		SELECT rid, (ST_DumpAsPolygons(raster)).geom AS geometrie, (ST_DumpAsPolygons(raster)).val
		FROM srtm.srtm_lr_wgs84_block
			JOIN dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text, 
					'SELECT code_dept, st_transform((geometrie), 4326) as geometrie
					FROM ign_bd_topo.departements_lr 
					WHERE code_dept=''34'''
					) dept_34(code_dept character varying, geometrie geometry)
			ON ST_Intersects(raster, geometrie)
		) srtm_34
	WHERE val>=600
	) AS zone_montagne_34
	JOIN dblink('hostaddr=host port=5432 dbname=db user=use password=pass'::text, '
				SELECT id_entite_occsol, area, perimeter, st_transform(cleangeometry(occsol_2006.geometrie), 4326) as geometrie,
					id_poste
				FROM occupation_du_sol.occsol_2006
					JOIN ign_bd_topo.departements_lr ON st_intersects(occsol_2006.geometrie, departements_lr.geometrie)
				WHERE id_poste LIKE ''31%'' AND code_dept=''34'''
				) ocsol_foret_34(id_entite_occsol integer, area double precision, 
					perimeter double precision, geometrie geometry, 
					id_poste character varying)
		ON ST_Intersects(zone_montagne_34.geometrie, ocsol_foret_34.geometrie)

SUCCES!! Temps d'exécution : près de 350 minutes, soit près de 6H !!! Possibilité d'améliorer l'exécution? A voir.

Pour info/rappel, les sous requêtes s'exécutent respectivement en 5 et 7.5 minutes

Cette requête a été retentée en modifiant la sous-requête ocsol_foret_34 en retirant l'intersection avec le département 34, en pensant que la dblink allait pouvoir utiliser les index de la table ocsol.

Résultat : temps d'exécution 380 minutes. Aucun gain!

Vectorisation du MNT du département de l'Hérault pour les altitudes >= à 600 m

Avec un MNT en tuiles de 250 px
SELECT val, geometrie
	FROM	(
		SELECT rid, (ST_DumpAsPolygons(raster)).geom AS geometrie, (ST_DumpAsPolygons(raster)).val
		FROM bd_topo.mnt_34_tif_block
		) mnt_34
	WHERE val>=600

temps d'éxécution : 39 minutes

Avec un MNT en tuiles de 60px
SELECT val, geometrie
	FROM	(
		SELECT rid, (ST_DumpAsPolygons(raster)).geom AS geometrie, (ST_DumpAsPolygons(raster)).val
		FROM bd_topo.mnt_34_tif_block60
		) mnt_34
	WHERE val>=600

temps d'exécution : 38 minutes

La taille des tuiles ne semble pas être un élément d'optimisation des requêtes.

Création d'une fonction pour récupérer les valeurs d'un raster d'une base distante à partir d'une couche vecteur

Récupération pour des données ponctuelles des altitudes correspondantes du MNT du département de l'Hérault stocké sur un autre serveur

Création d'une fonction opérant dans la base distante qui comporte les différents types et fonctions “raster”. Le paramètre de la fonction est une géométrie au format texte. 

CREATE OR REPLACE FUNCTION cen_calcul_mnt_point(text)
  RETURNS double precision AS
$BODY$	
DECLARE
	var_geom alias FOR $1;
	val_rec RECORD;
	var_alt double precision;
BEGIN
	FOR val_rec IN
	SELECT val
	FROM dblink('hostaddr=host port=5432 dbname=db user=use password=pass',
			'SELECT (gv).val as val
			FROM (
				SELECT ST_Intersection(raster, ST_Geomfromtext('''||var_geom||''', 2154)) AS gv
				FROM bd_topo.mnt_34_tif_block60
				WHERE ST_Intersects(raster, ST_Geomfromtext('''||var_geom||''', 2154))
			) foo'
		)mnt_test (val double precision)
	LOOP
	var_alt:= val_rec.val;
	END LOOP;
RETURN var_alt;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;

Exemple d'utilisation sur un échantillon de 10 données. 

SELECT id_obs, nom_complet, cen_calcul_mnt_point(ST_AsText(st_transform(geometrie, 2154))) AS alt
FROM saisie.saisie_observation
LIMIT 10

Résultat 

id_obs;nom_complet                                        ;alt
2     ;"Zerynthia polyxena (Denis & Schiffermüller, 1775)";
3     ;"Pyrgus malvae (Linnaeus, 1758)"                   ;59
28    ;"Pararge aegeria (Linnaeus, 1758)"                 ;64
17235 ;"Issoria lathonia (Linnaeus, 1758)"                ;
21248 ;"Chazara briseis (Linnaeus, 1764)"                 ;88
21690 ;"Gonepteryx cleopatra (Linnaeus, 1767)"            ;
21733 ;"Papilio sibilla Linnaeus, 1767"                   ;388
22122 ;"Maculinea arion Linnaeus, 1758"                   ;
22123 ;"Maculinea alcon rebeli Hirschke, 1904"            ;
22148 ;"Maculinea alcon D., 1775"                         ;

temps d'exécution : 1.6 s

NB: Si pour certaines données vecteurs nous n'avons pas de donnée d'altitude, c'est que ces données ne concernent pas le département de l'Hérault

Récupération pour des données surfaciques des altitudes correspondantes du MNT du département de l'Hérault stocké sur un autre serveur

Moyenne, Maximum, Minimum et Écart type

Création d'un TYPE de donnée composite qui contiendra les données calculées du MNT concernant un polygone : identifiant du polygone (pour les calcul d'agrégation), altitude moyenne, altitude maximum, altitude minimum et écart type. 

CREATE TYPE enum_stat_mnt AS
   (id integer,
    moy numeric,
    max numeric,
    min numeric,
    ecart_type numeric);

Création d'une fonction opérant dans la base distante qui comporte les différents types et fonctions “raster”.

Les paramètres de la fonction sont l'identifiant du polygone et la géométrie au format texte. 

CREATE OR REPLACE FUNCTION cen_calcul_stat_mnt_polygone(integer, text)
  RETURNS enum_stat_mnt AS
$BODY$	
DECLARE
	var_id alias FOR $1;
	var_geom alias FOR $2;
	myrec enum_stat_mnt;
BEGIN
	FOR myrec IN
 
		SELECT var_id AS id, round(avg(val)) AS moy, max(val) AS max, min(val) AS min, 
			round(stddev(val)::numeric, 2) AS ecart_type
		FROM dblink('hostaddr=host port=5432 dbname=db user=use password=pass',
			'SELECT (gv).val as val				
			FROM (
				SELECT ST_Intersection(raster, ST_Geomfromtext('''||var_geom||''', 2154)) AS gv
				FROM bd_topo.mnt_34_tif_block60
				WHERE ST_Intersects(raster, ST_Geomfromtext('''||var_geom||''', 2154))
				) foo'
		) mnt_polyg_test (val numeric)
		GROUP BY var_id
	LOOP
	END LOOP;
RETURN myrec;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
 

Exemple d'utilisation sur une zone tampon de 10 m pratiquée sur l'échantillon précédent. 

SELECT id_obs, nom_complet, (stat).moy, (stat).max, (stat).min, (stat).ecart_type
FROM (
	SELECT id_obs, nom_complet, 
		cen_calcul_stat_mnt_polygone(id_obs, ST_AsText(st_buffer(st_transform(geometrie, 2154), 10))) AS stat
	FROM saisie.saisie_observation
	LIMIT 10
	)polyg_test

Résultat: 

id_obs;nom_complet                                        ;moy;max;min;ecart_type
2     ;"Zerynthia polyxena (Denis & Schiffermüller, 1775)";
3     ;"Pyrgus malvae (Linnaeus, 1758)"                   ;58 ;59 ;57 ;1.41
28    ;"Pararge aegeria (Linnaeus, 1758)"                 ;64 ;64 ;64 ;
17235 ;"Issoria lathonia (Linnaeus, 1758)"                ;
21248 ;"Chazara briseis (Linnaeus, 1764)"                 ;92 ;101;84 ;7.68
21690 ;"Gonepteryx cleopatra (Linnaeus, 1767)"            ;
21733 ;"Papilio sibilla Linnaeus, 1767"                   ;388;388;388;
22122 ;"Maculinea arion Linnaeus, 1758"                   ;
22123 ;"Maculinea alcon rebeli Hirschke, 1904"            ;
22148 ;"Maculinea alcon D., 1775"                         ;

temps d'exécution : 1.7 s

Médianes du jeu de données (maximum et minimum) et Médiane calculée

Création d'un TYPE de données composite qui contiendra les données calculées du MNT concernant un polygone : médiane maximum, médiane minimum, médiane calculée. 

CREATE TYPE enum_median_mnt AS
   (min_med numeric,
    max_med numeric,
    calc_med numeric);

Création d'une fonction opérant dans la base distante qui comporte les différents types et fonctions “raster”.

Mais avant ça, une fonction intermédiaire est nécessaire pour calculer les occurrences cumulées des différentes valeurs les paramètres de cette “sous-fonction” sont: la valeur en question, l'identifiant et la géométrie de l'objet. 

CREATE OR REPLACE FUNCTION cen_mnt_occurence_cumul(numeric, integer, text)
  RETURNS numeric AS
$BODY$	
DECLARE
	var_val alias FOR $1;	
	var_id alias FOR $2;
	var_geom alias FOR $3;
	myrec record;
BEGIN
	FOR myrec IN
		SELECT sum(occ) AS occ_cum 
		FROM (
			SELECT var_id AS id, val AS alt, count(val) AS occ
			FROM dblink('hostaddr=host port=5432 dbname=db user=use password=pass',
					'SELECT (gv).val as val			
					FROM (
						SELECT ST_Intersection(raster, ST_Geomfromtext('''||var_geom||''', 2154)) AS gv
						FROM bd_topo.mnt_34_tif_block60
						WHERE ST_Intersects(raster, ST_Geomfromtext('''||var_geom||''', 2154))
						) foo
					ORDER by val'
				) mnt_polyg_test (val numeric)
			GROUP BY id, val
			ORDER BY id, val
		) mnt_polygone_occ 
		WHERE alt <= var_val
	LOOP
	END LOOP;
RETURN myrec.occ_cum;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;

Puis la “fonction principale”, dont les paramètre sont l'identifiant du polygone et la géométrie au format texte: 

CREATE OR REPLACE FUNCTION cen_calcul_mediane_mnt_polygone(integer, text)
  RETURNS enum_median_mnt AS
$BODY$	
DECLARE
	var_id alias FOR $1;
	var_geom alias FOR $2;
	myrec enum_median_mnt;
BEGIN
	FOR myrec IN
		WITH mnt_polygone_occ AS (
			SELECT var_id AS id, val AS alt, count(val) AS occ
			FROM dblink('hostaddr=host port=5432 dbname=db user=use password=pass',
					'SELECT (gv).val as val			
					FROM (
						SELECT ST_Intersection(raster, ST_Geomfromtext('''||var_geom||''', 2154)) AS gv
						FROM bd_topo.mnt_34_tif_block60
						WHERE ST_Intersects(raster, ST_Geomfromtext('''||var_geom||''', 2154))
						) foo
					ORDER by val'
				) mnt_polyg_test (val numeric)
			GROUP BY id, val
			ORDER BY id, val
		)
 
		SELECT min(mnt_polygone_occ.alt) AS min_median, max(mnt_polygone_occ.alt) AS max_median, 
			round((min(mnt_polygone_occ.alt)+max(mnt_polygone_occ.alt))/2, 2) AS calc_median
		FROM mnt_polygone_occ
		WHERE cen_mnt_occurence_cumul(alt, var_id, var_geom)>=(SELECT sum(mnt_polygone_occ.occ)/2 AS n_demi FROM mnt_polygone_occ)
 
 
	LOOP
	END LOOP;
RETURN myrec;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
 

Exemple d'utilisation sur une zone tampon de 10 m pratiquée sur l'échantillon précédent. 

SELECT id_obs, nom_complet, (med).min_med, (med).max_med, (med).calc_med
FROM (
	SELECT id_obs, nom_complet, 
		cen_calcul_mediane_mnt_polygone(id_obs, ST_AsText(st_buffer(st_transform(geometrie, 2154), 10))) AS med
	FROM saisie.saisie_observation
	LIMIT 10
        )polyg_med_test

Résultat: 

id_obs;nom_complet                                        ;min_med;max_med;calc_med
2     ;"Zerynthia polyxena (Denis & Schiffermüller, 1775)";
3     ;"Pyrgus malvae (Linnaeus, 1758)"                   ;57     ;59     ;58.00
28    ;"Pararge aegeria (Linnaeus, 1758)"                 ;64     ;64     ;64.00
17235 ;"Issoria lathonia (Linnaeus, 1758)"                ;
21248 ;"Chazara briseis (Linnaeus, 1764)"                 ;88     ;101    ;94.50
21690 ;"Gonepteryx cleopatra (Linnaeus, 1767)"            ;
21733 ;"Papilio sibilla Linnaeus, 1767"                   ;388    ;388    ;388.00
22122 ;"Maculinea arion Linnaeus, 1758"                   ;
22123 ;"Maculinea alcon rebeli Hirschke, 1904"            ;
22148 ;"Maculinea alcon D., 1775"                         ;

temps d'exécution : 4.3 s

Création d'une fonction pour récupérer les valeurs de plusieurs rasters d'une base distante à partir d'une couche vecteur

Récupération pour des données ponctuelles des altitudes correspondantes au MNT du département concerné stocké sur un autre serveur

Création d'une fonction opérant dans la base distante qui comporte les différents types et fonctions “raster”. Le paramètre de la fonction est une géométrie au format texte et le numéro du département extrait du code_insee de la commune où est localisée la donnée.

Chaque table RASTER étant nommée avec son numéro du département, on concatène dans le dblink() le paramètre “numéro de département” dans le nom de la table appelée.

Seuls les MNT des départements de la région Languedoc-Roussillon ont été intégrés, d'où le WHERE. 

CREATE OR REPLACE FUNCTION cen_calcul_mnt_point_v2(text, text)
  RETURNS double precision AS
$BODY$	
DECLARE
	var_geom alias FOR $1;
	val_rec RECORD;
	var_alt double precision;
	var_dept alias FOR $2;
BEGIN
	FOR val_rec IN	
	SELECT val
	FROM dblink('hostaddr=host port=5432 dbname=db user=user password=pw',
			'SELECT (gv).val as val
			FROM (
				SELECT ST_Intersection(raster, ST_Geomfromtext('''||var_geom||''', 2154)) AS gv
				FROM bd_topo.mnt_'||var_dept||'_asc_block60
				WHERE ST_Intersects(raster, ST_Geomfromtext('''||var_geom||''', 2154))
			) foo'
		)mnt_test (val double precision)
	WHERE var_dept IN ('11', '30', '34', '48', '66')
	LOOP
	var_alt:= val_rec.val;
	END LOOP;
RETURN var_alt;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;

Exemple d'utilisation sur les données dont le champ “elevation” n'est pas renseigné. 

SELECT id_obs, nom_complet, substring(code_insee FROM 1 FOR 2),
	cen_calcul_mnt_point_v2(ST_AsText(st_transform(geometrie, 2154)), substring(code_insee FROM 1 FOR 2)) AS alt
FROM saisie.saisie_observation
WHERE saisie_observation.elevation IS NULL

temps d'exécution pour 139 données : 40 s

Test du tutoriel de Pierre racine

http://trac.osgeo.org/postgis/wiki/WKTRasterTutorial01 SRTM/données odonates-papillons

chargement du SRTM dans postgis

chargement du raster en tuiles carrées de 60px de côté. 

raster2pgsql.py -r G:\TEMP\SRTM\srtm_37_04.tif -t srtm.srtm_lr_wgs84_block -s 4326 -c -f raster -F -I -M -k 60x60| psql -d raster_lr -h localhost -U yoann

La table produite comprend 10201 enregistrements/tuiles. Propriétés/métadonnées de la table 

SELECT (md).*, (bmd).* 
 FROM (SELECT ST_Metadata(raster) AS md, 
              ST_BandMetadata(raster) AS bmd 
       FROM srtm.srtm_lr_wgs84_block LIMIT 1
      ) foo;
visualisation de l'étendue du raster dans QGIS

avec l'extension Fast SQL layer

Bien nommer les champs de la requête avec des alias. Saisir les alias des champs “identifiant” et “géométrique” dans la barre de dessous de la fenêtre Query.

Attention! QuantumGIS, contrairement à OpenJUMP qui est utilisé dans le tutoriel, nécessite un champ unique!

étendue de chaque tuile
SELECT rid, raster::geometry AS geometrie
FROM srtm.srtm_lr_wgs84_block

étendue de l'ensemble
SELECT 1 AS id, ST_Buffer(ST_Union(raster::geometry), 0.000001) AS geometrie
FROM srtm.srtm_lr_wgs84_block

NB: Comme le souligne le tutoriel, la fonction ST_Buffer() garantit la fusion des polygones d'emprise de chaque tuile qui n'auraient été que largement partielle (polygones non exactement jointifs!?)

vectorisation d'une tuile
SELECT (ST_DumpAsPolygons(raster)).geom, 
        (ST_DumpAsPolygons(raster)).val
 FROM srtm.srtm_lr_wgs84_block
 WHERE rid=3299
 ORDER BY val

Temps d'exécution : 890 ms

Un problème se pose avec QGIS et la fonction Fast SQL layer. Alors que les géométries résultantes sont justes, les valeurs attribuées sont erronées (expl: 1077739520 pour une valeur attendue de 29.

Croisement SRTM/region biogéographique/Corine Land Cover

Ce cas consiste à récupérer les habitats ouverts de la zone méditerranéenne situés en dessous de 700 mètres d'altitudes.

Intégration de la mosaïque de dalles

Contrairement à ce qui a été fait dans les cas précédents, les dalles SRTM concernées (au nombre de 4) sont intégrées en tuiles de 50*50. L'intégration par tuiles de 60*60 crée des bordures vides. 

création de la table et intégration de la première dalle
raster2pgsql.py -r G:\TEMP\SRTM\srtm_37_04\srtm_37_04.asc -t srtm.srtm_block_50 -s 4326 -c -f raster -F -I -M -k 50x50| psql -d dbname -h host -U user
insertion des autres dalles
raster2pgsql.py -r G:\TEMP\SRTM\srtm_38_04\srtm_38_04.asc -t srtm.srtm_block_50 -s 4326 -a -f raster -F -k 50x50| psql -d dbname -h host -U user
raster2pgsql.py -r G:\TEMP\SRTM\srtm_38_03\srtm_38_03.asc -t srtm.srtm_block_50 -s 4326 -a -f raster -F -k 50x50| psql -d dbname -h host -U user
raster2pgsql.py -r G:\TEMP\SRTM\srtm_37_03\srtm_37_03.asc -t srtm.srtm_block_50 -s 4326 -a -f raster -F -k 50x50| psql -d dbname -h host -U user

Etape 1: Création de la couche vectorielle de la zone méditerranéenne sous 700 mètre d'altitude

Le traitement sur l'ensemble de la couche region_biogeo semble être trop lourd. La requête ci-dessous plante systématiquement. 

SELECT '<700m'as alti, ST_Buffer(ST_Union(geometrie), 0.000001) AS geometrie
FROM (
	SELECT val, geometrie AS geometrie
	FROM	(
		SELECT rid, (ST_DumpAsPolygons(raster)).geom AS geometrie, (ST_DumpAsPolygons(raster)).val
		FROM srtm.srtm_block_50
			JOIN biogeographie.region_biogeo ON ST_Intersects(raster, ST_Transform(geometrie, 4326))
		WHERE gid=4
		) srtm_med
	
	WHERE val<700
	) srtm_med_700

Pour contourner ce problème, une couche de la zone biogéographique désirée (“méditérranée”) est créée et subdivisée en 4 partie.

La requête, appliquée à une seule sous-partie, comme ce qui suit 

SELECT '<700m'as alti, ST_Buffer(ST_Union(geometrie), 0.000001) AS geometrie
	FROM (
		SELECT val, geometrie AS geometrie
		FROM	(
			SELECT rid, (ST_DumpAsPolygons(raster)).geom AS geometrie, (ST_DumpAsPolygons(raster)).val
			FROM srtm.srtm_block_50
				JOIN biogeographie.region_biogeo_med_part ON ST_Intersects(raster, ST_Transform(geometrie, 4326))
			WHERE dom_part='4-1'
			) srtm_med	
		WHERE val<700
		) srtm_med_700_1

renvoie son résultat avec succès en 42 min.

Cette dernière est intégré dans un table PostGIS pour la suite

ATTENTION! Pour une raison non déterminée, cette même requête crash et éteint le serveur Postgres pour seulement l'un des 4 secteurs. En divisant en 2 ce dernier, le problème persiste sur l'une des moitié, mais avec cette fois ci une erreur GEOS de type UnionCascaded. Il aurait été possible de continuer à diviser la partie à problème jusqu'à isolement et identification du problème, mais au vu de grand nombre de lignes, ce travail aurait pu être fastidieux. Il a été décidé d'utiliser plutôt les outils de PostGIS.

En décomposant la requête, il semble que le problème vienne du ST_Union. Les autres fonctions d'union (ST_MemUnion et ST_UnaryUnion) ont été testé mais sans succès. La solution est apporté par ST_Collect.

Par contre le buffer sur ST_Collect fait planter le serveur

Etape 2: Découpage de la couche d'occupation du sol selon le périmètre obtenu précédemment.