Suite

Intersection de plusieurs couches dans PostGIS

Intersection de plusieurs couches dans PostGIS


J'ai besoin d'écrire une procédure stockée PL/pgSQL pour le calcul de l'intersection de plusieurs couches. Par exemple, pour trois couches A, B, C, cette fonction doit renvoyer la géométrie ABC :

La fonction prend en entrée l'identifiant des couches qui doivent être intersectées. Quelqu'un peut-il me donner des conseils pour construire cette fonction? Mes réflexions à ce sujet, je dois croiser la première couche et la deuxième, puis le résultat de cette intersection croise la troisième couche, etc.


S'il vous plaît donnez-moi des commentaires sur ma solution proposée:

CRÉER OU REMPLACER LA FONCTION fp_intersect(lids varchar) RETOURNE l'entier AS $$ DECLARE lid_new entier ; lndx entier := 1; lids_arr entier[]; BEGIN IF regexp_replace(lids, E's+',", 'g') ~ E'^-?d+$' THEN RETURN -1; END IF; SELECT nextval ('g_layer_lid_seq') INTO lid_new; lids_arr : = string_to_array(regexp_replace(couvercles, E's+',", 'g'), ','); SUPPRIMER LA TABLE SI EXISTE tmp_intersect; CREATE TEMPORARY TABLE tmp_intersect AS SELECT geom FROM g_lgeom WHERE lid = lids_arr[1]; WHILE lndx < array_length(lids_arr, 1) LOOP DROP TABLE IF EXISTS tmp; CREATE TEMPORARY TABLE tmp AS SELECT ST_Intersection(geom, g_next) AS geom FROM tmp_intersect JOIN ( SELECT geom AS g_next FROM g_lgeom WHERE lid = lids_arr[lndx+1] ) AS _ ON ST_Intersects(geom, g_next); lndx := lndx+1; SUPPRIMER LA TABLE SI EXISTE tmp_intersect; CREATE TEMPORARY TABLE tmp_intersect AS SELECT geom FROM tmp; FIN DE BOUCLE ; INSÉRER DANS g_lgeom(lid, geom) SELECTIONNER lid_new, (_.p_geom).geom DE (SELECT ST_Dump(geom) AS p_geom FROM tmp_intersect) AS _; RETOUR couvercle_nouveau; FIN $$ LANGUE plpgsql;

Que se passe-t-il si vous vérifiez les intersections de la zone de chevauchement avant le calcul de l'intersection et de la zone, cela peut réduire le temps.

Vous obtiendrez probablement plus de réponses à ce sujet sur gis.stackexchange.com.

Il y a plusieurs choses que vous pouvez faire.

Vous devez vous assurer d'obtenir ce premier filtrage des polygones qui se croisent réellement à l'aide de l'index.

Mettez un index gist sur la table avec de nombreuses géométries et utilisez st_dwithin(geom,500) au lieu de st_intersects sur les géométries mises en mémoire tampon. C'est parce que les géométries tamponnées ne peuvent pas utiliser l'index calculé sur les géométries non tamponnées.

De plus, vous dites que vous avez plusieurs polygones. S'il y a en fait plus d'un polygone dans chaque multipolygone, vous obtiendrez peut-être beaucoup plus de vitesse si vous divisez d'abord les polygones en polygones simples avant de créer l'index. Cela permettra à the.index de faire une bien plus grande partie du travail.

Il existe en fait une fonction dans postgis pour diviser même des polygones simples en plus petits morceaux pour la même raison. ST_Sous-diviser

Utilisez donc d'abord ST_Dump pour obtenir des polygones uniques :

Enfin la requête, quelque chose comme

Si cela est toujours lent, vous pouvez commencer à jouer avec ST_SubDivide.

Encore une chose. Si le multipolygone unique de la table b contient de nombreuses géométries, divisez-les également et placez-y également un index.

Cela pourrait être lent aussi après toutes ces choses. Cela dépend du nombre de points de sommet dans les polygones divisés qui se coupent réellement (et pour st_dwithin également du nombre de points de sommet dans les polygones avec des cadres de délimitation qui se chevauchent)

Mais maintenant, vous n'avez aucun index pour vous aider, cela devrait donc le rendre beaucoup plus rapide.


PyQGIS Problème d'ajout de certaines couches vectorielles à QGIS à partir de postgres

J'ai créé un outil pour trouver des couches de notre base de données PostGIS en fonction d'une correspondance de chaîne et les ajouter au projet QGIS actuel. Dans l'exemple ci-dessous, je recherche ‘landscape’ dans le nom de l'ensemble de données.

Problème : il ajoute certains ensembles de données au projet en cours, mais pas tous (travail sur 2).

Le code de base qui ajoute les couches dans: ## Je nourris tous les paramètres avec ‘P’ à la fin

Notez qu'aucune de ces couches n'avait à l'origine le champ ID au début et même alors, les mêmes 2 sur 4 fonctionnaient. Ensuite, j'ai déposé et créé (dans FME à l'aide du rédacteur PostGIS) pour tous les ensembles de données et en plaçant la colonne ‘id’ au premier plan (voir les paramètres ci-dessous) et les deux mêmes ensembles de données ne voulaient toujours pas être chargés.

Lors de l'ajout de la couche à partir de la fenêtre du navigateur QGIS, les connexions PostGIS les deux ensembles de données problématiques sont ajoutés avec succès, fonctionne également si j'utilise le gestionnaire de base de données.

Les 4 ensembles de données ont les mêmes types de données d'attribut/colonne et les mêmes noms de colonne, id, geom, autocad_layer, etc. Je ne peux donc voir aucune différence dans les données pour expliquer le problème. Tous les ensembles de données ont vu leur géométrie vérifiée et corrigée avant d'entrer dans la base de données. Ce sont de vraies tables, -pas des vues- contrairement à un post précédent.

Q1. Avez-vous une idée de la raison pour laquelle ces 2 ensembles de données ne se chargent pas dans l'outil python ?

Q2. Quelle est la différence entre l'interface graphique QGIS (panneau de navigateur/gestionnaire de base de données) ajoute le code de couche et l'équivalent python ?

2 réponses

Le problème est probablement le type de géométrie. PostGIS peut avoir des tables de géométrie génériques, fme peut les écrire mais QGIS ne peut pas les gérer sans que vous spécifiiez le type de géométrie.

Vérifiez le type de géométrie de la colonne spatiale. Assurez-vous que la case "Créer des colonnes spatiales génériques" n'est pas cochée, puis vous pouvez choisir un point / une polyligne / un polygone, etc.

Répondu il y a 5 mois par nielsgerrits avec 4 votes positifs

Pour toute autre lecture, j'ai trouvé d'autres circonstances ou des détails supplémentaires lorsque les couches postgis ne se chargeaient pas dans QGIS :

Lorsque les couches contiennent des géométries mixtes dans un seul jeu de données ou des géométries complexes telles que des collections ou des surfaces polyédriques. Vous pouvez voir plusieurs couches du même nom dans le panneau du navigateur QGIS, une pour le point, la ligne, le polygone. Vous pouvez utiliser Postgis SQL pour forcer les collections en arcs/courbes CAO multi-polygones ( ST_CollectionExtract) en lignes (ST_CurveToLine).

Lorsque des erreurs de géométrie étaient présentes comme des auto-intersections et que vous pouvez les corriger à l'aide de postgis SQL ST_MakeValid.

Lorsqu'il manquait un système de coordonnées et que le SRID est signalé comme 0. Pour nous, cela était dû à l'utilisation d'un CRS personnalisé qui n'a pas d'EPSG/SRID enregistré. Vous pouvez utiliser Postgis SQL UpdateGeometrySRID à condition que vous ayez ajouté votre CRS personnalisé à la table spatial_ref et que vous connaissiez le numéro qui s'y rapporte.

Si vous n'avez pas de chance, les données sources telles que la CAO présentent de nombreux problèmes de ce type et vous devrez peut-être nettoyer les données dans un ordre particulier. Ce qui a fonctionné pour moi, c'est :

  1. collection à multi polygone,
  2. courbe à ligne
  3. St rendre valide
  4. Mettre à jour le SRID de la géométrie (s'applique à la couche)
  5. Ajoutez un champ 'id' int en série comme premier champ d'attribut (si ce n'est déjà fait) et définissez-le comme clé primaire car QGIS n'aime pas ne pas l'avoir. Vous constaterez peut-être qu'il se charge, mais vous ne pourrez pas modifier le calque.

Gardez à l'esprit que St makevalid peut réparer les choses et replacer les éléments géographiques dans une collection, vous devrez donc peut-être exécuter à nouveau une ou plusieurs étapes.

J'ai implémenté ces étapes dans l'atelier FME via un appelant python pour la publication sur Postgis. Mais également créé des outils python QGIS pour nettoyer les données déjà présentes dans la base de données.


Combiner des données avec des couches de polygones qui se chevauchent partiellement ? (QGIS)

J'ai un ensemble de données avec les quartiers de Baltimore, puis une série d'ensembles de données pour les secteurs de recensement de Baltimore remontant aux années 1930. Mon objectif est de voir comment la population et le nombre d'unités ont varié au fil des ans pour chaque quartier. Cependant, les quartiers et les secteurs de recensement se chevauchent de manière imprécise, de sorte que les situations suivantes se produisent :

  1. La plupart des quartiers contiennent plusieurs secteurs de recensement.
  2. Certains secteurs de recensement chevauchent plusieurs quartiers.
  3. Certains secteurs de recensement contiennent plus d'un quartier complet.

Je pense que la meilleure façon de gérer cela sera d'agréger les chiffres de population par le rapport pondéré du confinement d'un secteur de recensement donné dans un quartier donné. Comme quelques exemples (voir l'image aussi)
Dans le recensement de 1940, Reservoir Hill (quartier) contient 100 % de CT 0132 (pop 5 070), 80 % de CT 0131 (pop 4 795), 20 % de CT 0133 (pop 6263), et un petit éclat (disons 2%) de CT 0134 (pop 3 934). Donc 5070 + (4795*.8) + (6263*.2) + (3934*.02)=10.237.

Parkview/Woodbrook est un quartier entier contenu dans le CT 0134 en 1940. Étant donné que peut-être 10 % du 0134 est contenu dans le quartier, pop 3 934 * 0,1 = 393.

Je sais que cela introduira un biais d'agrégation, mais je pense que c'est probablement beaucoup plus précis que de simplement trouver le point central de chaque secteur de recensement, puis d'affecter la population totale à tout quartier contenant le point central.

Existe-t-il un moyen raisonnablement faisable d'effectuer ce type d'analyse et de calculs avec QGIS ?

Voici en quelque sorte ce que j'envisage comme sortie de quasi-tableaux dans le fichier de formes "Neighborhood" :

Quartier 1940s_CT 1940s_CT_Pleine Pop 1940s_CT_Area_Percentage 1940s_CT_Estimated_Pop
Colline du réservoir G240510000131,G240510000132,G240510000133,G240510000134 4795,5070,6263,3934 .8,1,.2,.02 10237
Vue sur le parc G240510000134 3934 .1 393.4

C'est le problème le plus lié que j'ai trouvé jusqu'à présent, bien que je n'aie jamais utilisé PostGIS.

2 réponses

QGIS Couche virtuelle (SpatiaLite SQL):

Personnellement, j'aime laisser les données de base intactes et avoir Couche virtuelles en place avec des analyses, plutôt que des données dupliquées. Aussi, personnellement, je préfère SQL presque partout :

N'oubliez pas de remplacer les valeurs enveloppées dans <> par les identifiants de couche/champ réels.

Répondu il y a 2 mois par géozelot avec 1 vote positif

Expressions QGIS :

Utilisez ces quatre expressions pour créer de nouveaux champs (virtuels ou non) dans votre Neighborhood Shapefile :

N'oubliez pas de remplacer les valeurs enveloppées dans <> par les identifiants de couche/champ réels.


Temple PSM dans le SIG

But:
Trouvez les meilleurs polygones correspondants dans une série temporelle à l'aide d'une base de données spatiale PostGIS
Arrière-plan:
Pour le projet précédent, une comparaison des données sur le logement des années 1930 aux données actuelles et historiques du recensement américain était requise. La partie géographique des données du recensement a été obtenue à partir du National Historical Geographic Information System (NHGIS) qui produit des fichiers de formes téléchargeables avec un champ GISJOIN permettant une liaison pratique du fichier de formes aux tableaux de données de recensement également disponibles via NHGIS. Les données sur le logement des années 1930 se présentaient sous la forme d'une image d'une carte de Philadelphie de 1937 Home Owners' Loan Corporation (HOLC) disponible ici. Très largement, les quartiers classés par HOLC en fonction de la race/de la classe et ces notes ont ensuite été utilisées pour déterminer la probabilité d'une aide au prêt au logement. Il y a un débat si le processus HOLC a normalisé les pratiques de prêt à domicile basées sur la race, puis conduit à d'autres problèmes de logement tels que la redlining. Plus d'informations.

Traiter:
Les données cartographiques HOLC ont été extraites d'un SIG et un fichier de formes a été créé pour le classement du quartier (ce processus est souvent répertorié comme « numérisation en tête-à-tête » dans les offres d'emploi) voici un didacticiel pratique pour ArcMap et un autre utilisant QGIS. Shapefile et HOLC jpeg disponibles sur github.

Le schéma de classement des quartiers HOLC ne coïncide avec aucun arrangement de recensement actuel (ou historique), de sorte que les capacités d'interrogation spatiale de pgAdmin et du moteur de base de données PostgreSQL ont été utilisées pour faire correspondre les données numériques HOLC nouvellement créées avec les données géographiques du recensement décennal NHGIS.
Créez une base de données spatiale PostGIS et utilisez pgShapeLoader pour charger les fichiers de formes (données spatiales HOLC et NHGIS)

Après de nombreux essais et erreurs, le code suivant a fonctionné le mieux :

créer un chevauchement de table2010census comme
avec limité comme(
sélectionnez distinct us2010.*
de nous2010, holc1937
O ST_DWithin(us2010.geom, holc1937.geom, 0,25)
commander par us2010.gid
)
/* la section ci-dessus limite la requête des polygones du recensement de 2010 (ensemble des États-Unis) à ceux proches des polygones HOLC – l'entrée peut varier en fonction de la projection des fichiers*/

SELECTIONNER DISTINCT ON (limited.gid)
limited.gid comme us2010a_gid,
holc1937.gid comme holc_gid, holc1937.name, holc1937.grade,
st_area((st_transform(limited.geom,4326)::geography)) comme area_check,
ST_Area(ST_Intersection((st_transform(limited.geom,4326)::geography),
(st_transform(holc1937.geom,4326)::geography))) comme intersect_area_2010_meters,
(ST_Area(ST_Intersection((st_transform(limited.geom,4326)::geography),
(st_transform(holc1937.geom,4326)::geography)))/
st_area((st_transform(limited.geom,4326)::geography)))*100 comme percent_overlap,

/*calcule l'aire des polygones de recensement et HOLC retournés et le pourcentage de chevauchement */

limité.*
DE limité, holc1937
où limited.geom && ST_Expand(holc1937.geom,1) et ST_Area(ST_Intersection(limited.geom, holc1937.geom)) >0

/* rejoint les fichiers où les geom se chevauchent*/

ORDER BY limited.gid, ST_Area(ST_Intersection((st_transform(limited.geom,4326)::geography),
(st_transform(holc1937.geom,4326)::geography))) DESC

/* trie par la plus grande quantité de chevauchement, et en conjonction avec l'appel DISTINCT ON
renvoie uniquement la première instance du polygone de recensement qui chevauche le polygone cible */

Bien que cette méthode nécessite que l'utilisateur mette à jour et répète ce code pour chaque comparaison, les tables créées peuvent ensuite être exportées vers votre SIG ou QGIS peut se connecter directement à votre base de données si vous préférez cette option. Comme la requête a également créé un identifiant unique pour chaque année et polygone des données du recensement, un tableau de concordance/clé a été créé à l'aide de ces données, ce qui a permis de joindre les données du recensement (race, statut de logement, etc.) le champ NHGIS GISJOIN mentionné précédemment. Sur la base des données de recensement pluriannuelles jointes via les tableaux de concordance, des comparaisons des régions géographiques HOLC pourraient alors être effectuées indépendamment des frontières changeantes du recensement.


PostSIG

La méthode QGIS convient si vous n'avez pas beaucoup de polygones à calculer, mais si vous avez un grand nombre d'entités, le processus prendra soit beaucoup de temps, soit plantera (d'ailleurs ArcGIS ne serait pas différent).

PostGIS à la rescousse. Pour cette approche, vous créez d'abord une base de données spatiale et activez l'extension PostGIS avec la commande CREATE EXTENSION postgis. Ensuite, vous pouvez charger les fichiers de formes dans PostGIS à l'aide du chargeur de fichiers de formes fourni avec PostGIS, ou vous pouvez utiliser le gestionnaire de bases de données QGIS pour les charger. Pendant le processus d'importation, vous devez spécifier que les couches sont en NAD 83 en spécifiant le code EPSG correct, en changeant le SRID de 0 à 4269.

PostGIS n'a pas beaucoup de définitions de système de coordonnées projetées mondiales ou continentales, nous devrons donc en ajouter une pour l'Amérique du Nord Albers Equal Area à sa table de référence spatiale. Une visite rapide de Spatial Reference et une recherche de ce système donnent la définition, et nous pouvons obtenir une instruction PostGIS Insert que nous pouvons copier et coller dans une fenêtre de requête SQL dans notre base de données. Avant de l'exécuter, nous devons changer le numéro SRID dans la déclaration de 9102008 à 102008 pour éviter de violer une restriction de contrôle qui empêche les ID d'être plus grands que 6 chiffres.

Une fois la définition en place, nous créons une série de tables vierges qui contiendront nos deux couches, puis exécutons une instruction insert dans laquelle nous prenons les colonnes que nous voulons des tables d'origine et les intégrons dans les nouvelles tables. Ce faisant, nous transformons également la géométrie de NAD 83 en Albers. À la fin, il est important de créer un index spatial sur la géométrie, car cela accélérera vraiment les sélections spatiales.

Une fois les données insérées, nous pouvons vérifier la validité de la géométrie avec ST_IsValid, et s'il y a une mauvaise géométrie, nous pouvons la corriger avec une autre instruction en utilisant ST_MakeValid, où IN contient des identifiants pour la mauvaise géométrie découverte dans l'instruction précédente.

Nous pouvons exécuter l'opération de chevauchement avec une seule instruction. PostGIS vous permet de calculer la surface à la volée avec la fonction ST_Area, et il existe deux fonctions pour le chevauchement : ST_Intersects agit comme une jointure spatiale qui relie une couche à l'autre en sélectionnant toutes les entités qui se croisent, tandis que ST_Intersection sélectionne les pièces réelles de chaque la géométrie des caractéristiques qui se coupent. Cet exemple concerne uniquement la Pennsylvanie, que nous sélectionnons à l'aide du code FIPS d'état 󈧮’ de la couche CD. C'est une bonne idée d'obtenir la bonne instruction sur un échantillon d'enregistrements avant de l'exécuter sur l'ensemble de l'ensemble. Les doubles deux-points sont un raccourci PostgreSQL pour convertir les types de données d'un type à l'autre. Ceci est nécessaire lorsque vous utilisez la fonction ROUND pour produire un résultat non entier (car ROUND ne peut pas être utilisé pour arrondir des nombres décimaux réels produits à partir de la fonction AREA à un nombre fixe de décimales).

Cette déclaration m'a pris environ 20 secondes à exécuter. Les résultats (voir ci-dessous) incluent plusieurs enregistrements que QGIS n'a pas renvoyés, où la zone et le chevauchement sont de 0, soit en raison d'une zone de chevauchement infiniment petite qui s'arrondit à zéro ou d'une interprétation stricte de l'intersection (qui inclut les zones qui se chevauchent et se touchent ). Bien qu'il existe une fonction ST_Overlap, elle ne retournera pas les géométries où une géométrie est complètement contenue dans une autre (nous ne pouvons donc pas l'utiliser). Par exemple, les ZCTA 19138 et 19139 apparaissent dans un district mais il y a deux enregistrements pour eux, un avec une valeur de 100 % et un autre avec une valeur de 0 %.

Résultat des opérations d'intersection et des calculs de surface dans pgAdmin / PostGIS

Nous pouvons jeter ces enregistrements en les supprimant du résultat final lorsque le processus est terminé, ou nous pouvons ajouter une autre instruction à notre clause WHERE pour les filtrer :

Cela a allongé le temps d'exécution à 30 secondes et fait passer le nombre d'enregistrements de 2 523 à 2 061.

Une fois que la déclaration semble correcte, nous pouvons supprimer le filtre AND pour la Pennsylvanie et générer un résultat pour l'ensemble du pays. En utilisant pgAdmin 4, nous pouvons écrire le résultat directement au format CSV. Ou, vous pouvez faire précéder l'instruction avec CREATE VIEW Chevauchement AS pour enregistrer l'instruction en tant que requête que vous pouvez appeler à tout moment. Ou, vous pouvez faire précéder l'instruction avec CREATE TABLE chevauchement AS et le résultat de la requête sera enregistré dans une nouvelle table. Cela prend plus de temps que les deux autres options, mais vous donne la possibilité d'interroger et de modifier la table résultante. L'exportation du tableau au format CSV peut être effectuée rapidement, vous offrant la meilleure des options 1 et 3. Le code final et le résultat sont indiqués ci-dessous.

Résultat final dans PostGIS / pgAdmin


Plus de messages

PostGIS 3.1.2

L'équipe PostGIS est heureuse de publier la version de PostGIS 3.1.2 !

Cette version est une version de correction de bogues, résolvant les problèmes trouvés dans la version 3.1 précédente.

    , TopoGeometry::geometry cast renvoie NULL pour les objets TopoGeometry vides (Sandro Santilli) , postgistigregéocodeur Meilleures réponses lorsqu'aucun zip n'est fourni (Regina Obe), gérer des dystems de coordonnées composées plus complexes (Paul Ramsey), ne faire que des retournements d'axe sur CRS qui ont un “Lat” comme première colonne (Paul Ramsey)
  • Prend en charge les versions récentes de Proj qui ont supprimé pjobtenirrelease (Paul Ramsey) , Ajuster la tolérance pour les calculs géodésiques (Paul Ramsey) , Conversion incorrecte de l'azimut géographique négatif en positif (Paul Ramsey) , Cluster DBSCAN non formé lorsque la longueur du jeu d'enregistrements est égale à minPoints (Dan Baston) , Mettre à jour les bboxes après scale/affine changements de coordonnées (Paul Ramsey) , Correction des problèmes de raster liés aux changements de tablefunc de PostgreSQL 14 (Paul Ramsey, Regina Obe) , mingw64 PostGIS / PostgreSQL 14 compile (Regina Obe, Tom Lane) , Mise à jour pour prendre en charge Tiger 2020 (Regina Obe) , Change Proj durée de vie du cache pour durer aussi longtemps que la connexion (Paul Ramsey), ajouter le support de build Pg14 (Paul Ramsey)

PostGIS 3.1.1

L'équipe PostGIS est heureuse de publier la version de PostGIS 3.1.1 !

Cette version est une version de correction de bogues, résolvant les problèmes trouvés dans la version 3.1 précédente.

    , Crash lors du passage de la collection avec uniquement des composants vides à ST_MakeValid , Faire fonctionner le pilote synthétique VSICURL comme documenté , Résultats instables de ST_MakeValid , Éviter de répertorier la même géométrie dans différentes collections

PostGIS 3.1.0

L'équipe PostGIS est heureuse de publier la version de PostGIS 3.1.0 !

Cette version expose les nouvelles fonctionnalités de GEOS 3.9 ainsi que de nombreuses améliorations de performances de base pour les jointures spatiales, l'accès aux objets volumineux, la sortie au format texte et plus encore.

Les performances sont une caractéristique clé de cette version, avec des améliorations apportées aux jointures spatiales, aux sorties de texte, aux lectures d'objets volumineux, à la sortie de tuiles vectorielles et à une multitude d'ajustements plus petits.

Le code de clustering k-means a été amélioré pour prendre en charge la pondération et les clusters de dimension supérieure.

Des générateurs de géométrie pour créer des pavages hexagonaux et carrés ont été ajoutés, pour des requêtes de résumé plus simples dans la base de données.

Enfin, PostGIS expose les dernières améliorations de la bibliothèque de géométrie GEOS version 3.9. Le nouveau moteur de superposition (alias “OverlayNG”) offre une gestion plus robuste des géométries d'entrée difficiles, en utilisant un ensemble de nouvelles stratégies de nœud pour traiter la géométrie. Pour l'utilisateur final, cela devrait signifier plus d'"exceptions de topologie" lors de l'utilisation des fonctions d'union, de différence, d'intersection ou de différence symétrique. PostGIS expose également la nouvelle capacité de superposition à précision fixe via un paramètre de taille de grille supplémentaire sur ST_Intersection et les autres fonctions de superposition.


Here be dragons: PostGIS 2.0 ajoute la prise en charge de la 3D, des images raster et de la topologie

Le système de base de données géographiques PostGIS approche de sa version 2.0, une mise à jour majeure qui ajoute plusieurs nouvelles fonctionnalités importantes. PostGIS est une extension du gestionnaire de base de données PostgreSQL qui implémente des types de données et des fonctions spécifiques au SIG. Tout aussi important, cependant, est le fait qu'un large éventail d'autres projets de système d'information géographique (SIG) open source peuvent utiliser PostGIS comme source de données principale, y compris les applications et les serveurs GUI. Cette version 2.0 est sur le point d'être une étape importante non seulement pour des raisons de vitesse et de stabilité, mais parce qu'elle étend la base de données dans trois nouveaux domaines : les données raster, la topologie et la 3D.

PostGIS ajoute la prise en charge des primitives géométriques (points, lignes, polygones, ainsi que des "collections" et autres structures de données des trois), ainsi que la lecture, la transformation et l'écriture d'une vaste gamme de formats de données géospatiales standard. Les applications SIG nécessitent souvent des opérateurs spéciaux pour calculer les distances et les surfaces, les unions, les intersections et d'autres fonctions définies, ainsi que des types de recherches spécialisés.

PostGIS implémente cette fonctionnalité en adhérant à la norme Simple Feature Access for SQL de l'Open Geospatial Consortium, bien que le projet ne paie pas pour les tests de conformité nécessaires pour se présenter comme une implémentation officielle. PostGIS est en développement actif depuis 2001 et a constitué une bibliothèque remarquable d'applications SIG de support, y compris GRASS GIS, gvSIG et MapServer. Il existe même des connecteurs commerciaux disponibles pour lier des produits SIG propriétaires aux bases de données PostGIS.

Pourtant, pendant tout ce temps, PostGIS s'est concentré sur ce que vous pourriez appeler la géométrie 2D traditionnelle, basée sur des vecteurs. Il ne faut pas voir cela comme un faiblesse la majorité des applications SIG sont 2D et vectorielles. Les données vectorielles codent les entités géographiques sous forme de formes : points et polygones sur une projection cartographique, lignes reliant les entités, etc.

Raster

Il existe cependant de nombreuses données disponibles sous forme d'images raster, ainsi que de photographies aériennes et satellites, par exemple, qui peuvent être géocodées afin qu'elles puissent être correctement alignées et transformées pour s'adapter à une carte. Ces dernières années, la communauté PostGIS a travaillé sur la prise en charge des rasters via un module complémentaire nommé à l'origine WKT Raster, qui a ensuite été renommé PostGIS Raster.

Avec 2.0, sa fonctionnalité sera entièrement intégrée à l'application principale. Les images raster sont prises en charge dans des tables raster spéciales, qui peuvent être chargées à partir de n'importe quel format pris en charge par la bibliothèque d'abstraction de données géospatiales (GDAL) et exportées vers n'importe quel format pris en charge par GDAL. La liste des formats pris en charge ne cesse de s'allonger, mais pour le moment, le projet GDAL en recense plus de 120.

Naturellement, le chargement et l'exportation ne sont pas le vrai travail de prise en charge des images raster en plus des vecteurs, donc PostGIS 2.0 ajoute une multitude de fonctions pour analyser et opérer sur les données à l'intérieur des pixels. Les rasters peuvent être "extrudés" dans la géométrie (par exemple, en transformant des régions d'une couleur dans l'image raster en un polygone), moyennés et inspectés. Les rasters peuvent être exploités avec des fonctions existantes (comme le calcul de leur intersection avec des formes géométriques vectorielles). Il sera également possible d'éditer les rasters importés en place dans la base de données, ce qui ouvre la porte à toutes sortes de transformations.

La troisième dimension

Les données 3D ne manquent pas non plus qui pourraient s'avérer intéressantes à examiner, si les types de données pour les stocker et les fonctions pour les transformer sont disponibles. PostGIS 2.0 ajoute une prise en charge étendue de la 3D, en commençant par deux types de géométrie : les surfaces polygonales et les réseaux triangulaires irréguliers (TIN). Les surfaces polygonales sont exactement ce que vous imaginez : des surfaces tridimensionnelles définies par des polygones connectés. Les TIN définissent des surfaces entièrement avec des triangles, mais la taille des triangles est flexible et ressemble un peu à un maillage multi-résolution dans un programme de modélisation 3D. Un TIN peut utiliser très peu de triangles pour décrire des zones plus plates de la carte, et plus pour décrire des caractéristiques importantes. Dans les deux cas, la prise en charge de la nouvelle géométrie ne consiste pas seulement en les types de données de base, mais en aide les opérateurs à effectuer des tâches récurrentes telles que la recherche de zones (et de volumes) de régions dans le nouveau format.

En plus des nouveaux types de géométrie, les indices spatiaux existants ont été rendus compatibles 3D et une bibliothèque de fonctions 3D ajoutée. Cela permettra aux utilisateurs de SIG de calculer des distances en trois dimensions, de trouver des intersections 3D de lignes et de formes, et de renvoyer des cadres de délimitation 3D ou de calculer des éléments complexes comme les plus courts chemins 3D.

L'application la plus simple de la 3D dans les SIG consiste à modéliser les caractéristiques du paysage en trois dimensions, mais il y en a plus. De nombreux autres types de données géocodées pourraient être importés dans une base de données PostGIS 3D pour la modélisation et l'analyse. Envisagez de calculer la visibilité directe entre les bâtiments, la propagation des ondes radio ou les vecteurs de vol, par exemple. Tous impliquent des opérations SIG familières, en 3D, même s'il ne s'agit pas strictement de tâches liées à la « cartographie ».

La meilleure façon de produire des données 3D est une question qui reste en grande partie le domaine de l'application frontale (nous y reviendrons un peu plus tard. ), mais PostGIS 2.0 ajoute la prise en charge de la sortie directe de données 3D au format X3D basé sur XML. X3D est défini par le Web3D Consortium, qui travaille dur pour que le format soit accepté dans HTML5. Que HTML5 autorise ou non l'inclusion de scènes X3D en ligne dans le contenu de la page, cependant, le format est susceptible d'être pris en charge dans les éléments <canvas> ou en tant qu'objets intégrés.

Topologie

La dernière extension majeure de la fonctionnalité PostGIS dans la version 2.0 est la prise en charge de la topologie. En termes de SIG, la topologie est ne pas une référence aux cartes topographiques. Cela signifie plutôt prendre en charge des types de données et des fonctions qui implémentent des graphiques mathématiques. En d'autres termes, au lieu de points, de lignes et de polygones, la topologie utilise des nœuds, des arêtes et des faces &mdash potentiellement même des arêtes dirigées et/ou pondérées.

La transformation de la géométrie vectorielle en données topologiques transforme ce qui était une collection de formes en une représentation mathématique de la scène, y compris la façon dont les nœuds et les régions sont connectés. L'œil humain peut établir ces connexions instantanément, mais la base de données nécessite un support explicite pour les exprimer. L'intégration de la version topologique d'une couche cartographique dans la base de données permet à l'application d'effectuer une recherche de chemin, un routage et une analyse de flux qui ne peuvent pas être effectués sur une géométrie brute.

L'exemple canonique de cette distinction est le problème du pont de Königsberg : un mathématicien humain (disons, Euler, pour en prendre un au hasard. ) peut regarder un dessin de la ville et de ses sept ponts célèbres et élaborer son graphe sous forme de nœuds et d'arêtes implicitement. . Une application PostGIS ne peut pas faire la même chose avec juste le formes de la rivière et des masses continentales. Les données doivent d'abord être converties.

PostGIS 2.0 sera capable de transformer la géométrie standard en données topologiques, de valider la topologie et de modifier les nœuds et les arêtes. La topologie peut également être convertie en langage de balisage géographique (GML) standard pour la sortie.

Cette étape représente le début de la prise en charge de la topologie dans PostGIS. Les futures applications incluent les flux de réseau (comme la modélisation du trafic), la gestion de crise pour la planification des catastrophes (comme les arbres couvrants minimum et les chemins les plus courts pour la logistique), la coloration automatique des cartes, les problèmes de « couverture » ​​et bien plus encore. Ceux d'entre nous qui ne travaillent pas dans le domaine des SIG peuvent avoir tendance à penser au travail SIG à l'échelle de la cartographie nationale, par exemple, mais considérer à quel point le support de la topologie serait précieux pour PostGIS pour l'analyse des flux de réseau ou pour le routage des services publics dans un bâtiment.

Prise en charge des applications et améliorations supplémentaires

Les fonctionnalités telles que la topologie et la 3D ont une valeur limitée sans prise en charge dans les applications qui utilisent PostGIS. Sur ce front, GRASS et gvSIG sont les premiers à inclure la topologie, d'autres devraient suivre. MapServer et QGIS, en revanche, prennent déjà en charge la couche raster. Il peut y avoir plus dans la catégorie de prise en charge des données raster qui héritent de leur fonctionnalité raster de la prise en charge antérieure du plug-in PostGIS Raster, la documentation n'est pas toujours claire.

Pour les données 3D, la seule application qui semble être prête à être prise en charge à l'heure actuelle est gvSIG, qui devrait l'inclure dans la prochaine version. Il s'agira d'une prise en charge de la visualisation 3D uniquement. Cependant, plusieurs blogs SIG open source sont enthousiasmés par les possibilités de la 3D, y compris son potentiel d'intégration avec la réalité virtuelle ou la réalité augmentée. Cela semble cependant être à quelques pas, tout comme la possibilité de combiner la topologie et les données 3D.

La prise en charge du raster, de la 3D et de la topologie prévue pour PostGIS 2.0 ne sont pas les seules nouvelles fonctionnalités du système de base de données. Les blogs et les listes de diffusion mentionnent à plusieurs reprises les améliorations apportées au géocodeur TIGER, qui importe des données provenant des données cartographiques du domaine public recueillies par le US Census Bureau, et qui s'est cassé lorsque TIGER a changé de format en 2010. Il existe également une nouvelle fonction de « géocodeur inversé ». qui prend un point de la carte et renvoie les données d'adresse ou les intersections de rues à proximité, et un chargeur de fichiers de forme GUI remanié qui, pour la première fois, peut charger plusieurs fichiers à la fois. Enfin, il existe des versions expérimentales pour Windows pour la toute première fois &mdash dans les versions précédentes, les utilisateurs de Windows devaient compiler PostGIS à partir de zéro. Le projet a une liste complète des fonctions nouvelles, améliorées et obsolètes dans sa documentation en ligne.

Il pourrait encore y avoir d'autres améliorations à venir. PostGIS n'aurait pas encore déclaré son gel final des fonctionnalités, et il existe des sous-traitants financés travaillant sur certains domaines importants, y compris la 3D. La version finale est officiellement attendue dans (l'hémisphère nord) "au début du printemps", qui n'a bien sûr que quelques jours à ce stade. En attendant, les audacieux peuvent obtenir des packages de test sur le site PostGIS. Avec ces nombreuses nouvelles fonctionnalités substantielles, il peut être intéressant de jeter un coup d'œil tôt.

Entrées d'index pour cet article
Articles invitésWillis, Nathan

(Connectez-vous pour poster des commentaires)

Here be dragons: PostGIS 2.0 ajoute la prise en charge de la 3D, des images raster et de la topologie

Publié le 31 mars 2011 à 3:56 UTC (jeu.) par flewellyn (abonné, #5047) [Lien]


4 réponses 4

J'essaierais de créer un tampon négatif, s'il mange des polygones minces, alors c'est bon, s'il ne mange pas le polygone, alors c'est le mien. :-)

exécutez ce script, après avoir défini au préalable 2/3 de la largeur des polygones linéaires .

créer la table name_table en tant que
SELECT ST_Buffer(
(ST_Dump(
(ST_Union(
ST_Tampon(
(geom),-0.0001))))).geom,
0,0001)) en tant que geom de source_table

OS :-).

en fin de compte, votre suggestion est ce qui a fonctionné le mieux pour moi. I ended using something like ST_Area(ST_Buffer(geom, -10)) , the -10 being -10 meters in my case. If anything returned 0 from that expression then I could filter it out.

Thank you, I am glad that it helped you, use this approach in such cases, with respect .

Instead of area/perimeter, it is better to use the area divided by the square of the perimeter (or its inverse).

This is also called "shape index". The square of the perimeter divided by the area has a minimum value of 4*Pi() (in the case of a disk, which is the most compact 2D geometry), so it can be normalized by 4*Pi() for an easy interpretation (normalized values close to 1 then mean that you have very compact objects and squares have a values of approximately 1.27).

EDIT: A threshold on the area would be usefull to remove the very small artefacts, which could be compact. Then the shape index would show better contrast.
EDIT: in addition to this answer, the use of ST_Snap could help you solve the problem before it occurs.

Merci! But I'm unsure how ST_Snap could help in this case. If I got it right, you're suggesting something like (o.overlap_perimeter^2 / o.overlap_area) / (4 * Pi()) as overlap_ratio ? This is having worse results for me than just area / perimeter.

Now using o.overlap_perimeter / (4 * sqrt(o.overlap_area)) as overlap_ratio according to this paper, but still worse results (although that's hard to quantify what I mean by worse) isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/I-7/135/…, page 183.

Thank you for this, I had never heard of the "shape index". I had always thought that using a minimum bounding rectangle was the best way to answer this sort of question. I found this, repository.asu.edu/attachments/111230/content/…, which is interesting.

– John Powell
Mar 20 at 16:50

@JohnPowell intersting paper, thanks. I see that what I know as a shape index is called circularity index in the paper. My problem with minimum bounding rectangles is that it doesn't work with very concave objects (e.g. U-shaped)

@bplmp ST_Snap would help you snap the vertices of "nearly" adjacent polygons so that they do no overlap anymore. There is no scale on your figures, but your artefact look like lines, so I guess that you can use a tolerance value theat is enough to avoid artefacts but does not affect the large polygons.

One option would be to use the ratio of the area of the polygon to the longest line that can be drawn using its extremities. Identifying long narrow polygons.

select * from polygons where ST_Length(ST_LongestLine(geom, geom)) < ST_Area(geom) * 4

This works pretty well for sliver polygons. You can adjust what the ratio (what you multiply the area with) to suit your needs and projection.

It sounds like this might match your use case: Eliminate selected polygons


Combines selected polygons of the input layer with certain adjacent polygons by erasing their common boundary. The adjacent polygon can be either the one with the largest or smallest area or the one sharing the largest common boundary with the polygon to be eliminated.

Eliminate is normally used to get rid of sliver polygons, i.e. tiny polygons that are a result of polygon intersection processes where boundaries of the inputs are similar but not identical.


It sounds like you'd want to try the "Largest Common Boundary" option.

I realize now you were asking for postgis solutions not qgis solutions. My apologies, I don't think postgis has an equivalent function but I'll leave this up for posterity.


Chapter 7 The cartographer as mediator: Cartographic representation from shared geographic information *

The advent of the Spatial Data Infrastructure (SDI) and distributed, agent-based computing present new challenges for the development of cartographic systems. The ability to effectively share meaning between system elements (semantic interoperability) remains an area of active research. The authors suggest a concept of geographic mediation (geomediation) that may be useful in framing current and emerging processes of integrating and representing geographic information in the context of distributed information. It is proposed that this concept can be used to guide practice and theory related to the design and evaluation of cartographic systems that can be seen as central to concepts like Taylor's cybercartography. The chapter presents a high-level architecture, the “open cartographic framework”, based on a mediator approach. Implications for using such an approach are discussed in this chapter. The chapter concludes by suggesting that an effective geomediation process necessitates consideration of both formal and negotiated processes for establishing its meaning.

This chapter is based on a PhD comprehensive paper presented as part of the lead author's PhD program at Carleton University, Ottawa, Ontario, Canada.


Voir la vidéo: PostGIS Databases and Geoprocessing