Présentation du

Système d'Information

start - [Bienvenue !]

admin_wiki_sicen — Thu, 02 Feb 2012 17:39:11 +0100

Download Flash Player 9 or better for full experience.

Bienvenue !

Ce wiki est mis en place pour présenter le SIG développé au CEN L-R depuis fin 2006 sur la plate-forme libre PostgreSQL/PostGIS pour en améliorer la documentation, en favoriser la diffusion et l'utilisation et pourquoi pas son amélioration collective.

Cette base de données spatiales permet la gestion et l'exploitation de données relatives :

à la faune, la flore et aux habitats naturels
aux limites des collectivités
aux périmètres de protection et d'inventaire
aux observateurs et aux partenaires
aux études réalisées et aux protocoles employés
aux référentiels taxonomiques
aux référentiels habitats
aux référentiels cartographiques vecteurs

La particularité de ce SIG repose depuis 2006 sur l'utilisation d'un système de gestion de bases de données relationnelles (utilisation des contraintes de domaine, de colonne, de table) pour la gestion de données spatiales et à thématique environnementale. L'utilisation de SGBD spatiaux est aujourd'hui de plus en plus répandue mais mérite de se généraliser¹⁾.

La mise en œuvre logicielle de ce système d'information repose sur l'utilisation de logiciels open source. La base de données est gérée par le couple PostgreSQL/PostGIS.

Vous trouverez sur ce site :

une présentation générale du SIG du CEN (ici)
les attentes qui ont guidées la mise en oeuvre de ce SI (ici à venir)
les choix logiciels effectués (ici)
la description des différents modules de la base de données (ici)
les codes des fonctions plpgsql et plpython développées et utilisées pour exploitation de la base (ici)
les applications logicielles libres utilisées pour mettre en œuvre ce SI (ici)

Mathieu BOSSAERT (sig à cenlr point org) Responsable du Système d'Information du CEN L-R

— Mathieu Bossaert 2011/08/18 12:34

¹⁾ Système de Gestion de Bases de Données Relationnelles

integration_donnees:externes:ign - [Calcul de la taille sur disque des tables du schéma inpn ...

admin_wiki_sicen — Wed, 01 Feb 2012 18:56:27 +0100

Intégration des thèmes de l'IGN à notre base de données

octobre 2010

Calcul de la taille sur disque des tables du schéma ign_bd_topo (5 départements)

pg_total_relation_size(regclass) renvoit l'espace disque total utilisé par la table de cet OID ou nom, incluant les index et les données de la partie TOAST

voici les statistiques pour notre serveur (RHEL 4, diques en ext4) :

SELECT * FROM (
SELECT f_table_name, pg_total_relation_size('ign_bd_topo.'||f_table_name), pg_size_pretty(pg_total_relation_size('ign_bd_topo.'||f_table_name))
  FROM geometry_columns
  WHERE f_table_schema='ign_bd_topo'

Au total, l'ensemble des tables prend une place d'environ 2.9 Giga, le mnt n'y est pas encore intégré (http://trac.osgeo.org/postgis/wiki/WKTRaster|piste]])

nom de la table	pg_total_relation_size()	pg_size_pretty(pg_total_relation_size())
aire_triage	172032	168 kB
bati_indifferencie	265838592	254 MB
bati_industriel	15687680	15 MB
bati_remarquable	2703360	2640 kB
canalisation_eau	155648	152 kB
chef_lieu	385024	376 kB
chemin	185024512	176 MB
cimetiere	794624	776 kB
conduite	90112	88 kB
construction_lineaire	19144704	18 MB
construction_ponctuelle	1294336	1264 kB
construction_surfacique	466944	456 kB
emprise	417792	408 kB
gare	90112	88 kB
hydronyme	475136	464 kB
lieu_dit	34013184	32 MB
lieu_dit_habite	4743168	4632 kB
lieu_dit_non_habite	4349952	4248 kB
ligne_electrique	434176	424 kB
ligne_orographique	31105024	30 MB
oronyme	2719744	2656 kB
pai_administratif_militaire	548864	536 kB
pai_culture_loisirs	385024	376 kB
pai_espace_naturel	4685824	4576 kB
pai_gestion_eaux	311296	304 kB
pai_hydrographie	475136	464 kB
pai_industriel_commercial	958464	936 kB
pai_orographie	2719744	2656 kB
pai_religieux	507904	496 kB
pai_sante	122880	120 kB
pai_science_enseignement	393216	384 kB
pai_sport	90112	88 kB
pai_transport	638976	624 kB
pai_zone_habitation	5169152	5048 kB
piste_aerodrome	172032	168 kB
point_eau	1253376	1224 kB
poste_transformation	114688	112 kB
pylone	2088960	2040 kB
reservoir	1941504	1896 kB
reservoir_eau	1662976	1624 kB
route	509304832	486 MB
route_nommee	130162688	124 MB
route_primaire	12091392	12 MB
route_secondaire	270352384	258 MB
surface_activite	2908160	2840 kB
surface_eau	17031168	16 MB
surface_route	483328	472 kB
terrain_sport	1474560	1440 kB
toponyme_communication	237568	232 kB
toponyme_divers	622592	608 kB
toponyme_ferre	81920	80 kB
transport_cable	106496	104 kB
troncon_cours_eau	104742912	100 MB
troncon_laisse	892928	872 kB
troncon_voie_ferree	2121728	2072 kB
zone_vegetation	1361412096	1298 MB
total	3008372736	2869 MB

juillet 2010

Création de la table des lieu-dits pour l'interface de saisie

Les toponymes sont stockés dans 4 tables de la BD TOPO :

lieu_dit_habite
lieu_dit_non_habite
oronyme
hydronyme

De manière à pouvoir proposer efficacement une aide à la saisie d'un lieu-dit selon l'endroit ou la donnée a été localisée sur la carte, nous allons rassembler ces informations dans une seule et même table nommée lieu_dit, dans le schéma public de notre base de données :

lieu_dit.sql

CREATE TABLE lieu_dit
(
  id character varying(24) PRIMARY KEY,
  origin_nom character varying(13),
  nom character varying(70),
  importance character varying(2),
  nature character varying(24),
  code_insee character varying(5), 
  /* Ce dernier champ n'existe pas initialement dans les tables "nymes" de la BD TOPO 
  Nous le rajoutons pour permettre des jointures avec la table commune et accélérer 
  le temps de récupération des lieu-dits proches d'un objet créé sur la carte */
  gid serial UNIQUE -- utile pour QGIS
)
WITH (
  OIDS=FALSE
);
 
/* création de la colonne géométrique de type point en lambert 93 */
 
SELECT AddGeometryColumn('lieu_dit', 'geometrie', 310024140, 'POINT', 2 );
 
/* remplissage de la table */
 
INSERT INTO lieu_dit(id, origin_nom, nom, importance, nature, geometrie)
SELECT id, origin_nom, nom, importance, nature, geometrie FROM ign_bd_topo.lieu_dit_habite;
INSERT INTO lieu_dit(id, origin_nom, nom, importance, nature, geometrie)
SELECT id, origin_nom, nom, importance, nature, geometrie FROM ign_bd_topo.lieu_dit_non_habite;
/*INSERT INTO lieu_dit(id, origin_nom, nom, importance, nature, geometrie)
SELECT id, origin_nom, nom, importance, nature, geometrie FROM ign_bd_topo.hydronyme;
INSERT INTO lieu_dit(id, origin_nom, nom, importance, nature, geometrie)
SELECT id, origin_nom, nom, importance, nature, geometrie FROM ign_bd_topo.oronyme;*/
 
/* indexation des géométries */
 
CREATE INDEX lieu_dit_geometrie_gist
  ON lieu_dit
  USING gist
  (geometrie);
 
/* mise à jour du code insee */
 
UPDATE lieu_dit
SET code_insee=ign_bd_topo.commune.code_insee
FROM ign_bd_topo.commune
WHERE st_intersects(ign_bd_topo.commune.geometrie, lieu_dit.geometrie);
 
/* indexation du code_insee des lieu-dits */
 
CREATE INDEX ind_lieu_dit_code_insee
  ON lieu_dit
  USING btree
  (code_insee);

Comment lister les lieux-dits situés dans un rayon de 2 km de l'objet créé

SELECT ign_bd_topo.commune.nom||' | '||lieu_dit.nom||' - '||round(st_distance(centroid(st_transform(saisie_faune.geometrie,310024140)), lieu_dit.geometrie))||' m.', round(st_distance(centroid(st_transform(saisie_faune.geometrie,310024140)), lieu_dit.geometrie)), ign_bd_topo.commune.nom
FROM (SELECT * FROM saisie_faune WHERE gid=4) AS saisie_faune JOIN ign_bd_topo.commune ON st_intersects(ign_bd_topo.commune.geometrie, centroid(st_transform(saisie_faune.geometrie,310024140)))
JOIN lieu_dit USING (code_insee)
WHERE st_intersects(st_buffer(st_transform(saisie_faune.geometrie,310024140),2000), lieu_dit.geometrie)
ORDER BY round(st_distance(st_transform(saisie_faune.geometrie,310024140), lieu_dit.geometrie)) ASC, ign_bd_topo.commune.nom

SELECT ign_bd_topo.commune.nom||' | '||lieu_dit.nom||' - '||round(st_distance(centroid(st_transform(geomfromtext(),310024140)), lieu_dit.geometrie))||' m.', round(st_distance(centroid(st_transform(geomfromtext(),310024140)), lieu_dit.geometrie)), ign_bd_topo.commune.nom
FROM ign_bd_topo.commune JOIN lieu_dit USING (code_insee) 
WHERE st_intersects(ign_bd_topo.commune.geometrie, centroid(st_transform(geomfromtext(),310024140)))
AND st_intersects(st_buffer(st_transform(geomfromtext(),310024140),2000), lieu_dit.geometrie)
ORDER BY round(st_distance(st_transform(geomfromtext(),310024140), lieu_dit.geometrie)) ASC, ign_bd_topo.commune.nom

février 2010

integration de la BD PARCELLAIRE dans le schema ign_bd_parcellaire

format : shp
projection : lambert 93

janvier 2010

Intégration de la BDTOPO

Transformation des MIF en shp

La BD topo a été initialement commandée en MIF/MID…

#!/bin/bash
set -x
i=1 # on initialise le compteur
option=''
for rep in `ls -d */`; do
	for mif in `find $rep -iname "*.MIF"`; do
                #echo $i : $PWD : $rep : $(echo $(basename $mif) | cut -d"." -f1) | tr '[A-Z]' '[a-z]'
                echo ogr2ogr $option -a_srs \"+init=IGNF:LAMB93\" -f \"ESRI Shapefile\" SHP/$(echo $(basename $mif) | cut -d"." -f1 | tr '[A-Z]' '[a-z]').shp ./$mif
                #echo ogr2ogr $option -f \"PostgreSQL\" PG:\"host=localhost user=dba dbname=ign password=dbroot\" $PWD/$rep$mif -nln bdtopo.$(echo $mif | cut -d"." -f1) | tr '[A-Z]' '[a-z]' 
		eval ogr2ogr $option -a_srs \"+init=IGNF:LAMB93\" -f \"ESRI Shapefile\" SHP/$(echo $(basename $mif) | cut -d"." -f1 | tr '[A-Z]' '[a-z]').shp ./$mif
	done
	let $[ i+=1 ] # incremente i de 1 a chaque boucle
	option='-update -append' # après la première série de fichier, il faudra ajouter les nouvelles données aux fichiers shp déja existants.
done

Intégration des shp à la base de données

#!/bin/bash
for shp in `ls -R *.shp`; do
	eval shp2pgsql -d -W "latin1" -s 310024140 -g geometrie -I -S $shp ign_bd_topo.$(echo $(basename $shp) | cut -d"." -f1 | tr '[A-Z]' '[a-z]')>../SQL/$(echo $(basename $shp) | cut -d"." -f1 | tr '[A-Z]' '[a-z]').sql
	eval psql -f ../SQL/$(echo $(basename $shp) | cut -d"." -f1 | tr '[A-Z]' '[a-z]').sql -h localhost -U dba sicen 2>>../erreurs.log
done

février 2009

Transformer les shapes en sql : ”“for /R %F IN (*.shp) DO shp2pgsql -s 27572 -g geometrie -I ”%F” ”%~nF”>“ign_bdcarto_%~nF.sql””“

Chargement dans la base de données : ”“for /R %F IN (*.sql) DO psql -f ”%F” -h addresse_ip_hote -U utilisateur base_de_donnée”“

- Données de la BD Parcellaire (Département du Gard)

octobre 2008

intégration des différentes couches de la BD CARTHAGE

La version mise à disposition sur le site du sandre introduit des accents dans le nommage des colonnes.

Après téléchargement des couches (ici au format mif/mid, si shp alors passer à l'étape suivante), les transformer en shp à l'aide de la console (windows):

for /R %F IN (*.MIF) DO ogr2ogr -f "ESRI ShapeFile" -s_srs "EPSG:2154" -t_srs "EPSG:27582" "%~nF.shp" %F

Transformer ces shapes en sql :

for /R %F IN (*.shp) DO shp2pgsql -s 27582 -S -g geometrie -I  "%F"  "ign_bdc_%~nF">"prefixe_%~nF_suffixe.sql"

Les charger dans la base de données :

for /R %F IN (*.sql) DO psql -f "%F" -h addresse_ip_hote -U utilisateur base_de_donnée

cours:bts:session2 - [Saisie de données]

admin_wiki_sicen — Wed, 01 Feb 2012 14:16:27 +0100

BTS GPN - Systèmes d'information géographique

Session n°2 - SIG : structuration, création et consultation de données

Montpellier - le 1er février 2012

Mathieu Bossaert : sig@cenlr.org

Rappels de la première session

Projections
Cordonnées
Formats de données
Données attributaires
Opérateurs

Rappels sur l'utilisation de QGis

SIG : Création de données "vecteurs"

Structuration de la table attributaire
- Choix des descriteurs de la données (attributs)
Choix du type de géométrie

SIG : Saisie de données

pointage sur fond de carte
pointage de coordonnées
gps
aide à la saisie
- activation de l'édition topologiques (anti-recouvrement et options d'accrochage)

Exercice 1

Création d'un projet de cartographie sur les salins du castellas

Enregistrer le projet vide sous tp_session_2.qgs

Spécifier la projection du projet (lambert 93) et activer la projection à la volée

préférences → propriétés du projet →activer la projection à la volée → 2154→ trouver → ok

Ajouter au projet les couches raster (scan25 et orthophotos)

Précisez leur projection si les données ne s'affichent pas au bon endroit

* reprendre les deux étapes précédentes

Ajouter la couche ZPS du serveur WFS de la DREAL

Activer l'extension WFS (Web Feature Service)

Extensions → Gestionnaire d'extensions → cocher “Extension WFS” → Ok

Créer un nouveau serveur WFS

http://ws.carmen.developpement-durable.gouv.fr/WFS/19/dreal_lr_general

Se connecter au serveur et choisir la couche ZPS

Exercice 2

Création de couches de données vecteur en Lambert 93

Ajouter au fichier issu du GPS les colonnes nécessaires

ajouter le fichier gpx au projet
enregistrer la couche en shp
regarder les colonnes disponibles
nous allons supprimer les colonnes superflues en utilisant l'extension “table manager”
ajouter les colonnes manquantes (espece, abondance, unité)

Création de 3 couches de données

création d'une table de données d'observations de faune (points)
- id entier, date texte, espece texte, observateurs texte, abondance texte, unite texte, remarque texte
création d'une table de digitalisation des chemins (ligne)
- id numerique entier, nom texte 100, largeur entier
création d'une table de cartographie des bassins visibles (polygones)
- id numerique entier, nom texte 100

Création d'une nouvelle couche

On peut utiliser les menus ou la combinaison des touches Ctrl + Shift + n

On choisit le type de données et la projection

On ajoute les attributs (colonnes)

L'étape n° 3 est utile pour les attributs de type texte, pour préciser le nombre de caractères que le champ peut contenir. Répétez ces 3 ou 4 étapes pour chaque attribut. Quand la liste des attributs est complète, cliquez su “Ok” et données un nom à votre couche.

Si vous souhaitez par la suite modifier votre couche et ajouter un attribut par exemple, il existe l'extension Table Manager

Saisie de données

Choisir la couche à éditer (1) et basculer en mode édition (2)

Pour commencer à créer des données : cliquez sur capturer la ligne

On termine la saisie d'une nouvelle donnée par un clic sur le bouton droit de la souris. A ce moment, la fenêtre de renseignement des attributs s'ouvre.

Création d'une carte - Composeur d'impression

Installation du plugin OpenLayers

Il permet l'affichage d'un fond issue de google ou de yahoo

Installation

Aller dans extensions / installateur d’extensions python / Dans l'onglet “dépots” :

cliquez sur “ajouter un dépôt-tiers…”

Dans l'onglet “Extensions” :

Saisissez “openlayers plugin” dans la case filtre
sélectionnez le plugin qui apparait dans la fenêtre
puis cliquez sur installer l'extension

Une fois l'extension télécharger : dans les préférences du projet : Préférences / propriétés du projet /Système de coordonnées de référence /

Cocher la case activer la projection à la volée
Rechercher la valeur 3644
Cliqez sur Ok

Enfin, ajouter l'image de fonds en allant dans “Extension” / “Openlayers plugin” / “Add google Satellite layer”

Pour aller plus loin

Dans les propriétés du projet, vous pouvez spécifier un certain nombre de préférences qui vous permettront une saisie propre et facilitée :

définir les options d'accrochage permettra lors de la création d'un nouvel objet de voir votre curseur se coller à l'objet voisin
si vous créez par exemple une couche de polygones d'occupation du sol, il n'y a aucun raison que deux polygones voisins se chevauchent, vous pouvez choisir d'éviter les intersections entre nouveaux polygones

Accueil
Nuage de tags
derniers ajouts
Actualités
Présentation diverses
Offres de stage
Démarche
Création de la base
Les Données
- Diffusion / Echanges
  - Pistes de reflexion
- Format de données
  - Format Standard de données
- Intégration de données
  - Données internes
    - Problématique
  - Données externes
    - Généralités
    - DGI
    - DIREN L-R
    - IFN
    - IGN
    - INPN
    - SIGLR
- Sources de données
Exploitation de la base
- Commande d'export
- Contenu de la bdd
(Biblio|Webo)graphie
Outils utilisés
- SIG
  - PostgresSQL/PostGIS
    - Ressources
    - Fonds d'écran /Wallpapers
  - PostGIS/Raster
    - Ressources
    - Cas d'utilisation
  - Moskitt/Moskitt
    - Ressources
    - Cas d'utilisation
  - GRASS
  - QGis
  - GVSIG
  - gdal
    - Ressources
    - Cas d'utilisation
  - ogr2ogr
    - Ressources
    - Cas d'utilisation
  - Talend Open Studio / SDI
  - MDWEB2
  - Mapproxy
    - Installation
- Non SIG
  - iReport
  - Anuko Time Tracker
  - Expressions régulières
  - JabRef
    - Ressources
    - Cas d'utilisation
  - Dokuwiki
    - Ressources
    - Cas d'utilisation
  - Drupal
    - Cas d'utilisation
- Système d'exploitation
  - Linux
    - commandes Shell
  - RHEL / CentOS
  - FEDORA
    - Ressources en ligne
Cours
- BTS
  - Session 1
  - Session 2
bric à brac
- remplacement de codes html
- optimisation de pdf
FCEN
- CTT du 28 janvier 2011
- outil de suivi des actions de gestion
Interne
- gap_analysis
  - ressouirces
Confidentialité
modeles_de_page
- pages de fonction

outils:postgis_raster:cas_d_utilisation - [Etape 1: Création de la couche vectorielle de la ...

yoyo — Fri, 27 Jan 2012 10:03:53 +0100

PAGE EN COURS DE CONSTRUCTION

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

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)

non tuilée: 6001×6001 px
tuiles de 60×60 px : 10 201 tuiles

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

non tuilée: 5401×3481 px
tuiles de 250×250 px : 308 tuiles
tuiles de 60×60 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 .

 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;

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 :

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

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

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

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

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

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.

cours:bts:session1 - [Projection Lambert : conique]

admin_wiki_sicen — Wed, 25 Jan 2012 12:43:01 +0100

BTS GPN - Systèmes d'information géographique

Session n°1 - Le SIG, définition, objectifs, concepts et outils

Montpellier - le 25 janvier 2012

Mathieu Bossaert : sig@cenlr.org

Objectifs de ces séances

~~Faire de vous des spécialistes du SIG~~
Montrer l'importance des SIG nos métiers
Démystifier le SIG
Vous montrer que ce n'est pas si compliquer
Vous donner des éléments pour approfondir vos compétences dans ce domaine
- liens vers des documentation
- liens vers des données
Vous proposer un premier contact avec des outils que vous serez amenés à utiliser

SIG : définitions

Outil

Le plus souvent, on parle de SIG pour désigner un logiciel de cartographie. C'est une utilisation réductrice du terme. Mais c'est bien de cela que nous allons parler.

Les besoins

Diagnostic de site
Faire des préconisations de gestion
Suivi des travaux de gestion
Inventaire d'espèces, d'habitats
Inventaire des équipements
- martelières
- platelage
Inventaire des chemins
…

Concepts

Représenter sur un plan, un objet sphérique

Concepts - Projection

Projection de Mercator : cylindrique

L’Amérique du Sud semble plus petite que le Groenland, en réalité, elle est neuf fois plus grande :

17,8 millions de km² contre 2,1 millions¹⁾.

Concepts - Projection

Projection Lambert : conique

Les logiciels de SIG s’accommodent très bien des différentes projections et pratiquent la reprojection à la volée. Importance du paramétrage des préférences du logiciel

TP configuration de la projection dans QGIS

Installation de QGIS

Format - Types de données : vecteur

Format - Types de données : image

Format - Types de données : grille

L'espace est découpé en une grille. Chaque partie de la grille (pixel) contient une valeur.

Concepts - Cordonnées

Latitude (lat.)

La latitude exprime le positionnement nord-sud d'un point sur Terre, au nord ou au sud de l'équateur.
La latitude est une mesure angulaire s'étendant de 0° à l'équateur à 90° aux pôles.

Longitude (long.)

La longitude exprime le positionnement est-ouest d'un point sur Terre.
La longitude est donc une mesure angulaire sur 360° par rapport à un méridien de référence,de -180° (ouest) à +180° (est).
Le méridien de référence est le méridien de Greenwich.

Concepts - Cordonnées

Les coordonnées géographiques sont traditionnellement exprimées dans le Système sexagésimal :

DMS Degré:Minute:Seconde (49 ° 30 ′ 00 ″ - 123 ° 30 ′ 00 ″ )
DM Degré:Minute (49 ° 30.0 ' - 123 ° 30.0 ')
DD Degré décimal (49.5000 ° - 123.5000 °), généralement avec 4 décimales.

Le WGS 84 est le système géodésique associé au système GPS il s'est rapidement imposé comme la référence universelle pour la cartographie.

Mapinfo
Arcview
Geoconcept

libres/gratuits

GRASS
Quantum Gis
GvSIG
Ilwis

Quantum GIS

Quantum GIS est un logiciel SIG publié sous licence GPL.
Né en mai 2002
projets officiels de la fondation OSGeo
Fonctionne sous Linux, Mac OS X, ainsi que Windows
prend en charge de nombreux formats vectoriels et rasters
Liste sans cesse accrue de fonctionnalités via un puissant système d'extensions.

ubuntu : http://doc.ubuntu-fr.org/qgis
fedora : http://hub.qgis.org/projects/quantum-gis/wiki/DownloadFr#Fedora
RHEL, CentOS, Scientific Linux : http://elgis.argeo.org/

plugins

Convertisseur de couches OGR: pour simplifier l'export de couches dans différents formats.
TableManager: Pour modifier une table attributaire.
Le plugin “OpenLayers”

occupation du sol
- Corine Land Cover
- ocsol LR
IGN
- photos anciennes de l'IGN, certaines gratuites
- données et logiciels gratuits
Périmètres d'inventaire et périmètres réglementaires
- DREAL L-R de votre région
  - http://ws.carmen.developpement-durable.gouv.fr/WFS/19/dreal_lr_general
  - http://ws.carmen.developpement-durable.gouv.fr/WMS/19/dreal_lr_general
- MNHN (mise à jour moins récente qu'en DREAL)
Divers
- NASA
- openstreetmap
- http://www.naturalearthdata.com/downloads/
- données cliamtiques (grilles)

Afficher les raster (images) du dossier “données\rasters\SCAN1000”. Il s'agit du fond de carte de l'IGN
Ajouter à la carte la couche des départements français : données\vecteurs\GEOFLA
Zoomer sur l'emprise de la couche (départements)
Ajouter la trace GPX des stations de limonium : données\vecteurs\BIODIV
Zoomer sur l'emprise de la couche

Exercice n°2

utilisation de ces données

afficher seulement les département de la région “Languedoc-Roussillon”
renommez cette sélection “Languedoc-Roussillon”
afficher en plus les autres départements de France
renommez cette sélection couche “autres régions”
afficher les départements de la couche “languedoc-Roussillon” avec un fond vert et un contour pointillé noir
faites en sorte d'afficher les départements de la couche “autres régions” de manière à ce que les département d'une même région aient la même couleur (une couleur par région)

¹⁾ http://fr.wikipedia.org/wiki/Projection_de_Mercator

Présentation duSystème d'Information