LAPT (Linux And Perl Tips)

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