LAPT (Linux And Perl Tips)

Je suis en train d'insérer des "zillions" de points provenant d'un campagne de mesure lidar dans une base de données postgresql/postgis et le taux d'insertion/sec est "abassourdissant". Bon il est vrai qu'il a fallu un peu "adapter" les paramètres du fichier de configuration postgresql.conf pour ces opérations d'insertions massives et il est aussi vrai que l'utilisation de la bonne méthode d'insertion est déterminante. Mais tout de même pour un GNU/Linux Debian 6.0.5 Squeeze 64 bit tournant dans une machine virtuelle kvm sur un host physique Proxmox "ça décoiffe sec" les amis. Tout commence par la compilation de la dernière version de liblas  la version 1.7.0 en ce moment. Ensuite on rentre dans le vif du sujet je convertis mes fichiers las en csv avec l'utilitaire las2txt  : las2txt --parse MxyzcipnreadRGB  --precision 3 3 3 3 -i 1243-14-d_color_ortho_r1.las -o 1243-14-d_color_ortho_r1.csv Ensuite avec un  w

Prendre la geomA et la geomB et retourner : “ Dimensionally Extended 9 Intersection Model (DE-9IM)” selon section 4.3.6 de la doc postgis SELECT ST_Relate('SRID=21781;POINT(538494 152932)'::geometry,'SRID=21781;POINT(538494 152932)'::geometry) --"0FF-FFF-FF2" --deux points diffrents SELECT ST_Relate('SRID=21781;POINT(538494 152932)'::geometry,'SRID=21781;POINT(538492 152932)'::geometry) --"FF0-FFF-0F2" --deux lignes identiques SELECT ST_Relate('SRID=21781;LINESTRING(53 15, 55 16)'::geometry,'SRID=21781;LINESTRING(53 15, 55 16)'::geometry) --"1FF-F0F-FF2" --deux lignes partiellment identiques SELECT ST_Relate('SRID=21781;LINESTRING(53 15, 55 16)'::geometry,'SRID=21781;LINESTRING(53 15, 55 16, 58 18)'::geometry) --""1FF-00F-102" --deux lignes qui se coupe perpendiculaire SELECT ST_Relate('SRID=21781;LINESTRING(0 0, 10 0)'::geometry,'SRID=21781;LINESTRING(5 5, 5 -

CREATE OR REPLACE FUNCTION cncggetisosurface(geomdepart geometry,tolerance float,altitude double precision) RETURNS geometry AS $$ DECLARE idcarredepart integer:=0; carredepart geometry; altituderef double precision:=0; tempidgrid integer:=0; tempgeom geometry; tempaltitude double precision:=0; toutlageometrieavant geometry; toutlageometrie geometry; isendofstory boolean :=false; BEGIN IF GeometryType(geomdepart) = 'POINT' THEN FOR tempidgrid,tempaltitude,tempgeom IN (SELECT idgrid_elevation,mns_lidar, the_geom FROM l3d_grid1m_elevation WHERE ST_DWITHIN(the_geom,geomdepart, 0.01) ) LOOP IF tempaltitude > altituderef THEN idcarredepart := tempidgrid; altituderef := tempaltitude ; carredepart := tempgeom; END IF; END LOOP; --au depart etait le carre toutlageometrie := carredepart; ELSE --ensuite etait la geom carre,polygone etc... toutlageometrie := geomdepart; END IF; WHILE NOT isendofstory LOOP

CREATE SEQUENCE seq_police_isolines_id MINVALUE 1 START WITH 1 INCREMENT BY 1 CREATE TABLE police_isolines as (SELECT nextval('seq_police_isolines_id') as id, COUNT(*), ST_UNION(the_geom), FLOOR(mns_lidar) FROM l3d_grid1m_elevation WHERE ST_Intersects(the_geom, (SELECT geom as gb from bati_pol where No_eca = '11330')) --AND the_geom && (SELECT geom as gb from bati_pol where No_eca = '11330') GROUP BY FLOOR(mns_lidar) ORDER BY 3 ) ;

première étape charger un mns dans postgis:   ./raster2pgsql -I -C -r -t 50x50 -d -s 21781 /home/xao/postgis/MNS_zone_test.asc mns >/tmp/mns.sql puis on insère dans notre base de données: psql -U postgres -f /tmp/mns.sql gis01 enfin on peut créer une table temporaire avec une séquence CREATE SEQUENCE seq_bati_isolines_id MINVALUE 1 START WITH 1 INCREMENT BY 1 CACHE 10 CREATE TABLE bati_isolines as (SELECT nextval('seq_bati_isolines_id') as id, no_eca, (gv).geom, (gv).val as z FROM (SELECT b.no_eca , ST_Intersection(geom,rast) gv FROM bati_pol b, mns WHERE no_eca like '11330' AND st_intersects(geom,rast)) as t ) une fois ce petit détail reglé on peut allez voir nos isolignes dans QGIS avec le plugin "DBManager" avec la requête SELECT id,geom,z FROM bati_isolines

Cette fois on veut la liste des distances entre tous les points composants un polygone par exemple pour vérifier qu'il n'y a aucun point qui s'approche d'un autre à moins d'une certaine distance WITH line as (SELECT st_geometryn(st_boundary(geom),1) as g FROM bati_pol WHERE no_eca like '11330'), series as (SELECT generate_series(1,ST_NumPoints (g) ) as n FROM line), points1 as (SELECT st_pointn(g,n) as p1, n as n1 FROM line,series), points2 as (SELECT st_pointn(g,n) as p2,n as n2 FROM line,series) SELECT n1,n2,st_distance(p1,p2) FROM points1,points2 WHERE n1 != n2 ORDER BY 3 ASC

variante même requête en utilisant les Command Table Expression WITH line as (SELECT st_geometryn(st_boundary(geom),1) as g FROM bati_pol WHERE no_eca like '11330'), series as (SELECT generate_series(1,ST_NumPoints (g) ) as n FROM line) SELECT 'point '|| n,st_astext(st_pointn(g,n)), 'point '|| n+1,st_astext(st_pointn(g,n+1)), st_distance(st_pointn(g,n),st_pointn(g,n+1)) FROM line,series

Comment calculer la distance entre les points successifs d'un polygone SELECT ST_AsText(ST_PointN((st_geometryn(st_boundary(geom),1)), generate_series(1,ST_NumPoints (st_geometryn(st_boundary(geom),1))-1) )), ST_AsText(ST_PointN((st_geometryn(st_boundary(geom),1)), generate_series(1,ST_NumPoints (st_geometryn(st_boundary(geom),1))-1)+1 )), ST_Distance( ST_PointN((st_geometryn(st_boundary(geom),1)),s.i ) , ST_PointN((st_geometryn(st_boundary(geom),1)), generate_series(1,ST_NumPoints (st_geometryn(st_boundary(geom),1))-1)+1 ) ) FROM bati_pol , (SELECT generate_series(1,ST_NumPoints (st_geometryn(st_boundary(bati_pol.geom),1))-1) FROM bati_pol WHERE no_eca like '11330') as s(i) WHERE no_eca like '11330'

Pour remplacer toutes les fonctions de créations géométries avantageusement par une seule syntaxe générique on utilise un string WKT avec un cast: 'SRID=21781;POINT(538460 152947)'::geometry par exemple dans le contexte d'une requête : SELECT * FROM l3d_grid1m_elevation g WHERE st_intersects(g.the_geom,(SELECT st_buffer('SRID=21781;POINT(538460 152947)'::geometry,1))) AND g.the_geom && (SELECT st_buffer('SRID=21781;POINT(538460 152947)'::geometry,1)) Ainsi on a quasiment un remplacement unique pour toutes les fonctions du §8.3 Geometry constructors du manuel postgis

Je viens d'être confronté à un besoin d'estimation rapide de la densité de points par m2 d'un fihier "las" provenant d'un vol Lidar . Il fallait pouvoir répondre à la question combien de zones de 1m2 contiennnent des densités de points inférieur à un certain seuil. Wikipedia l'explique bien mieux mais en résumé pour les pressés : La télédétection par laser ou LIDAR «  li ght d etection a nd r anging  » est une technique qui permet en substance de "mesurer" depuis un avion ou un hélicoptère la position d'objets aux sol et d'obtenir un nuage de points x,y,z  qui donne un modèle numérique de surface ou MNS. En fonction du plan vol, de la météo des conditions topologiques en présence (grand bâtiments) on peut avoir des "zones d'ombres" ou on a peu ou pas de  données, il faut néanmoins pouvoir quantifier et visualiser ces zones. Pour les curieux je recommande l'excellente présentation de M. Marc Riedo responsable du SITN s