<box round red 100%>**//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.
</box>
{{tag>PostGIS-Raster SRTM BDTOPO classification}}

====== Installation de PostGIS-Raster/PostGIS 2.0 pour Windows XP======
===== Installation de PostGIS 2.0 =====
  - Télécharger les binaires compatible avec la version de PostreSQL (http://www.postgis.org/download/windows/experimental.php)
  - 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
  - Lancer //makepostgisdb.bat//
  - 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 :
  * //height// : hauteur de l'image
  * //scale x/y// : taille du pixel en x ou y
  * //skew x/y// : distance en x et y due à la rotation
  * //srid// : identifiant du référentiel spatial
  * //upperleft x/y// : coordonnée x et y du coin haut gauche de l'image
  * //width// :  largeur de l'image
==== 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'.
<code>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</code>

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

<code sql>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;</code>

Temps d'exécution : 46 minutes

====== CAS D'UTILISATIONS ======

==Caractéristiques des images utilisées:==
**SRTM (WGS84)**   
        * non tuilée: 6001x6001 px
        * tuiles  de 60x60 px : 10 201 tuiles

**MNT de département de l'Hérault(Lambert93)**   
        * non tuilée: 5401x3481 px  
        * tuiles  de 250x250 px : 308 tuiles 
        * tuiles  de 60x60 px : 5 369 tuiles

===== 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 .

<code> 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</code>

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
<code>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
		)</code>
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)  =====

<code sql>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;
</code>

ou 

<code sql>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;
</code>

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

Quant à 
<code sql>
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</code>
cela renvoie l'erreur
<code>ERREUR: Error creating GEOS Coordinate Sequence
État SQL :XX000
Contexte : PL/pgSQL function "_st_intersects" line 28 at affectation</code>
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 :

  * Attribution de la valeur d'altitude du **MNT** au **point d'observation de 10 individu** issus d'une autre base (via //dblink//).

<code sql>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</code>
Résultat en 3.5 s:
<code>
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
</code>


  * Même chose mais avec une **zone tampon de 100 m** autour des données ponctuelles pour attribuer à chaque espèce les différentes altitudes environnantes.
//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//
<code sql>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</code>
Temps d'éxécution : 3.5 s
<code> 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
...     |...                                    |...    |...  
</code>

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

<code sql>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</code>
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 ==

<code sql>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)</code>
**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==
<code sql>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</code>
temps d'éxécution : 39 minutes

==Avec un MNT en tuiles de 60px==
<code sql>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</code>
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.

<code sql>
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;</code>
  
Exemple d'utilisation sur un échantillon de 10 données.
<code sql>
SELECT id_obs, nom_complet, cen_calcul_mnt_point(ST_AsText(st_transform(geometrie, 2154))) as alt
FROM saisie.saisie_observation
LIMIT 10
</code>

Résultat
<code>
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"                         ;
</code>
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.
<code sql>
CREATE TYPE enum_stat_mnt AS
   (id integer,
    moy numeric,
    max numeric,
    min numeric,
    ecart_type numeric);
</code>

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.
<code sql>
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;
 </code>

Exemple d'utilisation sur une zone tampon de 10 m pratiquée sur l'échantillon précédent.
<code sql>
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
</code>

Résultat:
<code>
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"                         ;
</code>
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.
<code sql>
CREATE TYPE enum_median_mnt AS
   (min_med numeric,
    max_med numeric,
    calc_med numeric);
</code>

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.
<code sql>
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;
</code>

Puis la "fonction principale",
dont les paramètre sont l'identifiant du polygone et la géométrie au format texte:
<code sql>
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;
 </code>

Exemple d'utilisation sur une zone tampon de 10 m pratiquée sur l'échantillon précédent.
<code sql>
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
</code>

Résultat:
<code>
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"                         ;
</code>
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.


<code sql>
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;</code>
  
Exemple d'utilisation sur les données dont le champ "elevation" n'est pas renseigné.
<code sql>
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
</code>

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é.
<code>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
</code>
La table produite comprend 10201 enregistrements/tuiles.
Propriétés/métadonnées de la table
<code sql>SELECT (md).*, (bmd).* 
 FROM (SELECT ST_Metadata(raster) AS md, 
              ST_BandMetadata(raster) AS bmd 
       FROM srtm.srtm_lr_wgs84_block LIMIT 1
      ) foo;</code>
== 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==
<code sql>SELECT rid, raster::geometry as geometrie
FROM srtm.srtm_lr_wgs84_block</code>
{{:outils:postgis_raster:fast_query_qgis_raster.png?300|}}
==étendue de l'ensemble==
<code sql>SELECT 1 as id, ST_Buffer(ST_Union(raster::geometry), 0.000001) as geometrie
FROM srtm.srtm_lr_wgs84_block</code>
//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!?) //

{{:outils:postgis_raster:fast_query_qgis_raster_union.png?300|}}
==vectorisation d'une tuile==
<code sql>SELECT (ST_DumpAsPolygons(raster)).geom, 
        (ST_DumpAsPolygons(raster)).val
 FROM srtm.srtm_lr_wgs84_block
 WHERE rid=3299
 ORDER BY val</code>
Temps d'exécution : 890 ms
  * **visualisation dans QGIS**
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//.
  * **visualisation dans OpenJUMP** 
{{:outils:postgis_raster:sql_raster_oj_dap.png?300|}}

====== 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.
<code>
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
</code>


===== 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.
<code>
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
</code>

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
<code>
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
</code>
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.=====