Changes between Version 23 and Version 24 of UsersWikiSimplifyPreserveTopology


Ignore:
Timestamp:
Apr 28, 2012, 3:59:11 AM (10 years ago)
Author:
NicolasRibotOsgeo
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • UsersWikiSimplifyPreserveTopology

    v23 v24  
    174174
    175175{{{
    176 -- takes a table name, gid and geom column and simplify the layer based on given tolerance.
    177 -- returns a setof (gid int, geom geometry). where geom is the simplified geometry and gid is the tablename gid.
    178 create or replace function simplifyLayerPreserveTopology (tablename text, id_col text, geom_col text, tolerance float) returns setof record as $$
     176-- Simplify the given table of multipolygon with the given tolerance.
     177-- This function preserves the connection between polygons and try to avoid generating gaps between objects.
     178-- To identify objects after simplification, area comparison is performed, instead of PIP test, that may fail
     179-- with odd-shaped polygons. Area comparison may also failed on some cases
     180
     181-- Example: (table 'departement' is in the public schema)
     182-- select * from simplifyLayerPreserveTopology('', 'departement', 'gid', 'geom', 10000) as (gid int, geom geometry);
     183--
     184-- @param schename: text, the schema name of the table to simplify. set to null or empty string to use search_path-defined schemas
     185-- @param tablename: text, the name of the table to simplify
     186-- @param idcol: text, the name of a unique table identifier column. This is the gid returned by the function
     187-- @param tolerance: float, the simplify tolerance, in object's unit
     188-- @return a setof (gid, geom) where gid is the identifier of the multipolygon, geom is the simplified geometry
     189create or replace function simplifyLayerPreserveTopology (schemaname text, tablename text, idcol text, geom_col text, tolerance float)
     190returns setof record as $$
    179191    DECLARE
    180         tabname alias for $1;
    181         tid alias for $2;
    182         geo alias for $3;
    183         tol alias for $4;
     192        schname alias for $1;
     193        tabname alias for $2;
     194        tid alias for $3;
     195        geo alias for $4;
     196        tol alias for $5;
    184197        numpoints int:=0;
    185198        time text:='';
     199        fullname text := '';
    186200
    187201    BEGIN
    188 
    189         EXECUTE 'select sum(st_npoints(geom)), to_char(clock_timestamp(), ''MI:ss:MS'') from '||quote_ident(tabname) into numpoints, time;
     202        IF schname IS NULL OR length(schname) = 0 THEN
     203            fullname := quote_ident(tabname);
     204        ELSE
     205            fullname := quote_ident(schname)||'.'||quote_ident(tabname);
     206        END IF;
     207       
     208        raise notice 'fullname: %', fullname;
     209
     210        EXECUTE 'select sum(st_npoints('||quote_ident(geo)||')), to_char(clock_timestamp(), ''MI:ss:MS'') from '
     211            ||fullname into numpoints, time;
    190212        raise notice 'Num points in %: %. Time: %', tabname, numpoints, time;
    191213       
    192         EXECUTE 'create unlogged table poly as ('
    193                 ||'select '||quote_ident(tid)||', (st_dump('||quote_ident(geo)||')).* from '||quote_ident(tabname)||')';
     214        EXECUTE 'create unlogged table public.poly as ('
     215                ||'select '||quote_ident(tid)||', (st_dump('||quote_ident(geo)||')).* from '||fullname||')';
    194216
    195217        -- extract rings out of polygons
    196218        create unlogged table rings as
    197         select st_exteriorRing((st_dumpRings(geom)).geom) as g from poly;
     219        select st_exteriorRing((st_dumpRings(geom)).geom) as g from public.poly;
    198220       
    199221        select to_char(clock_timestamp(), 'MI:ss:MS') into time;
     
    234256        -- distinct clause to eliminate overlaping segments that may prevent polygon to be created,
    235257        -- then dump the collection of polygons into individual parts, in order to rebuild our layer.
    236         create temp table simplepolys as
     258        drop table if exists simplepolys;
     259        create  table simplepolys as
    237260        select (st_dump(st_polygonize(distinct g))).geom as g
    238261        from simplelines;
     
    252275        -- as input set is multipolygon, we first explode multipolygons into their polygons, to
    253276        -- be able to find islands and set them the right departement code.
    254         RETURN QUERY EXECUTE 'select '||quote_ident(tid)||', st_collect(geom) as geom '
     277        RETURN QUERY EXECUTE 'select '||quote_ident(tid)||', st_collect('||quote_ident(geo)||') as geom '
    255278            ||'from ('
    256279            --||'    select distinct on (d.'||quote_ident(tid)||') d.'||quote_ident(tid)||', s.g as geom '
    257280            ||'    select d.'||quote_ident(tid)||', s.g as geom '
    258             ||'   from '||quote_ident(tabname)||' d, simplepolys s '
     281            ||'   from '||fullname||' d, simplepolys s '
     282            --||'    where (st_intersects(d.'||quote_ident(geo)||', s.g) or st_contains(d.'||quote_ident(geo)||', s.g))'
    259283            ||'    where st_intersects(d.'||quote_ident(geo)||', s.g) '
    260284            ||'    and st_area(st_intersection(s.g, d.'||quote_ident(geo)||'))/st_area(s.g) > 0.5 '
     
    262286            ||'group by '||quote_ident(tid);
    263287           
    264         drop table simplepolys;
     288        --drop table simplepolys;
    265289       
    266290        RETURN;