<pre>Hello!<br><br>This patch add GEOS simplify functions to PostGIS<br>(tested only with PostGIS 1.1.6 and geos-3.0.0rc2 but I can build path to PostGIS 1.2.0)<br><br>With it, you can use large tolerance values without getting null result values...
<br>SQL : simplify_geos(geometry,float8) = GEOS DouglasPeuckerSimplifier::simplify(g1, tolerance)<br>SQL : simplifytopologypreserve_geos(geometry,float8) = GEOS TopologyPreservingSimplifier::simplify(g1, tolerance)<br><br>
Example : I use simplify_geos into online SVG Map application (demo here : <a href="http://siti.cmm.qc.ca/">http://siti.cmm.qc.ca/</a>)<br>PostGIS 1.1.6/Geos-3.0.0rc2/PostgreSQL 8.2<br><br>Could this patch be included by default into PostGIS?
<br> <br>Thanks a lots for PostGIS!!!<br><br>François <br><br>diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos.c<br>/new/postgis-1.1.6/lwgeom/lwgeom_geos.c<br>*** /postgis-1.1.6/lwgeom/lwgeom_geos.c        Mon Dec 12 09:40:50 2005<br>
--- /new/postgis-1.1.6/lwgeom/lwgeom_geos.c     Fri Nov 17 08:30:15 2006<br>***************<br>*** 66,71 ****<br>--- 66,75 ----<br>   Datum linemerge(PG_FUNCTION_ARGS);<br>  Datum LWGEOM_buildarea(PG_FUNCTION_ARGS);<br><br>+<br>
+ Datum simplify_geos(PG_FUNCTION_ARGS);<br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS);<br>+<br><br>   typedef  struct Geometry Geometry;<br><br>***************<br>*** 101,106 ****<br>--- 105,115 ----<br>   extern Geometry *GEOSGetCentroid(Geometry *g, int *failure);
<br>   extern bool GEOSHasZ(Geometry *g1);<br><br>+<br>+ extern Geometry *GEOSSimplify(Geometry *g1, double tolerance);<br>+ extern Geometry *GEOSTopologyPreserveSimplify(Geometry *g1, double tolerance);<br>+<br>+<br>   extern void GEOSSetSRID(Geometry *g, int SRID);
<br>   extern void GEOSdeleteChar(char *a);<br>   extern void GEOSdeleteGeometry(Geometry *a);<br>***************<br>*** 2896,2898 ****<br>--- 2905,3089 ----<br>   PG_RETURN_NULL();<br>  }<br><br>+<br>+ PG_FUNCTION_INFO_V1(simplify_geos);
<br>+ Datum simplify_geos(PG_FUNCTION_ARGS)<br>+ {<br>+       PG_LWGEOM *geom1;<br>+    double tolerance;<br>+    GEOSGeom g1,g2;<br>+      PG_LWGEOM *result;<br>+   int is3d;<br>+    int SRID;<br>+<br>+ #ifdef PROFILE<br>+       profstart(PROF_QRUN);
<br>+ #endif<br>+<br>+        geom1 = (PG_LWGEOM *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));<br>+         tolerance = PG_GETARG_FLOAT8(1);<br>+<br>+  is3d = ( TYPE_HASZ(geom1->type) );<br>+     SRID = pglwgeom_getSRID(geom1);<br>+         initGEOS(lwnotice, lwnotice);
<br>+<br>+ #ifdef PGIS_DEBUG<br>+     elog(NOTICE,"simplify starting");<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+   profstart(PROF_P2G1);<br>+ #endif<br>+      g1 = POSTGIS2GEOS(geom1);<br>+ #ifdef PROFILE<br>+  profstop(PROF_P2G1);
<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+   elog(NOTICE," constructed geometry - calling geos");<br>+       elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));<br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */
<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+      profstart(PROF_GRUN);<br>+ #endif<br>+      g2 = GEOSSimplify(g1,tolerance);<br>+ #ifdef PROFILE<br>+   profstop(PROF_GRUN);<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE,"simplify finished");
<br>+ #endif<br>+<br>+        if (g2 == NULL)<br>+      {<br>+            elog(ERROR,"GEOS Simplify() threw an error!");<br>+             GEOSGeom_destroy(g1);<br>+                PG_RETURN_NULL(); /* never get here */<br>+       }<br>+ #ifdef PGIS_DEBUG<br>
+       elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+  GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE<br>+     profstart(PROF_G2P);<br>+ #endif<br>+       result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br>
+       profstop(PROF_G2P);<br>+ #endif<br>+<br>+     if (result == NULL)<br>+  {<br>+            GEOSGeom_destroy(g1);<br>+                GEOSGeom_destroy(g2);<br>+                elog(ERROR,"GEOS Simplify() threw an error (result postgis<br>geometry formation)!");
<br>+             PG_RETURN_NULL(); /* never get here */<br>+       }<br>+    GEOSGeom_destroy(g1);<br>+        GEOSGeom_destroy(g2);<br>+<br>+     /* compressType(result); */<br>+<br>+ #ifdef PROFILE<br>+     profstop(PROF_QRUN);<br>+         profreport("geos",geom1, result);
<br>+ #endif<br>+   PG_FREE_IF_COPY(geom1, 0);<br>+<br>+        PG_RETURN_POINTER(result);<br>+ }<br>+<br>+ PG_FUNCTION_INFO_V1(simplifytopologypreserve_geos);<br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS)<br>+ {<br>
+       PG_LWGEOM *geom1;<br>+    double tolerance;<br>+    GEOSGeom g1,g2;<br>+      PG_LWGEOM *result;<br>+   int is3d;<br>+    int SRID;<br>+<br>+ #ifdef PROFILE<br>+       profstart(PROF_QRUN);<br>+ #endif<br>+<br>+   geom1 = (PG_LWGEOM *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
<br>+     tolerance = PG_GETARG_FLOAT8(1);<br>+<br>+  is3d = ( TYPE_HASZ(geom1->type) );<br>+     SRID = pglwgeom_getSRID(geom1);<br>+         initGEOS(lwnotice, lwnotice);<br>+<br>+ #ifdef PGIS_DEBUG<br>+        elog(NOTICE,"simplify starting");
<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+      profstart(PROF_P2G1);<br>+ #endif<br>+      g1 = POSTGIS2GEOS(geom1);<br>+ #ifdef PROFILE<br>+  profstop(PROF_P2G1);<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE," constructed geometry - calling geos");
<br>+     elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));<br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+      profstart(PROF_GRUN);<br>+ #endif<br>+      g2 = GEOSTopologyPreserveSimplify(g1,tolerance);
<br>+ #ifdef PROFILE<br>+   profstop(PROF_GRUN);<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE,"simplify finished");<br>+ #endif<br>+<br>+     if (g2 == NULL)<br>+      {<br>+            elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error!");
<br>+             GEOSGeom_destroy(g1);<br>+                PG_RETURN_NULL(); /* never get here */<br>+       }<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+  GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE
<br>+     profstart(PROF_G2P);<br>+ #endif<br>+       result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br>+   profstop(PROF_G2P);<br>+ #endif<br>+<br>+     if (result == NULL)<br>+  {<br>+            GEOSGeom_destroy(g1);<br>+                GEOSGeom_destroy(g2);
<br>+             elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error (result<br>postgis geometry formation)!");<br>+                PG_RETURN_NULL(); /* never get here */<br>+       }<br>+    GEOSGeom_destroy(g1);<br>+        GEOSGeom_destroy(g2);
<br>+<br>+  /* compressType(result); */<br>+<br>+ #ifdef PROFILE<br>+     profstop(PROF_QRUN);<br>+         profreport("geos",geom1, result);<br>+ #endif<br>+        PG_FREE_IF_COPY(geom1, 0);<br>+<br>+        PG_RETURN_POINTER(result);
<br>+ }<br>+<br>+<br>diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos_c.c<br>/new/postgis-1.1.6/lwgeom/lwgeom_geos_c.c<br>*** /postgis-1.1.6/lwgeom/lwgeom_geos_c.c     Wed Oct 25 03:23:00 2006<br>--- /new/postgis-1.1.6/lwgeom/lwgeom_geos_c.c Fri Nov 17 08:15:59 2006
<br>***************<br>*** 67,72 ****<br>--- 67,76 ----<br>  Datum LWGEOM_buildarea(PG_FUNCTION_ARGS);<br>   Datum linemerge(PG_FUNCTION_ARGS);<br><br>+<br>+ Datum simplify_geos(PG_FUNCTION_ARGS);<br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS);
<br>+<br><br>   LWGEOM *GEOS2LWGEOM(GEOSGeom geom, char want3d);<br>   PG_LWGEOM *GEOS2POSTGIS(GEOSGeom geom, char want3d);<br>***************<br>*** 2974,2976 ****<br>--- 2978,3162 ----<br>          PG_RETURN_POINTER(result);<br>
<br>  }<br>+<br>+<br>+ PG_FUNCTION_INFO_V1(simplify_geos);<br>+ Datum simplify_geos(PG_FUNCTION_ARGS)<br>+ {<br>+     PG_LWGEOM *geom1;<br>+    double tolerance;<br>+    GEOSGeom g1,g2;<br>+      PG_LWGEOM *result;<br>+   int is3d;<br>
+       int SRID;<br>+<br>+ #ifdef PROFILE<br>+       profstart(PROF_QRUN);<br>+ #endif<br>+<br>+   geom1 = (PG_LWGEOM *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));<br>+         tolerance = PG_GETARG_FLOAT8(1);<br>+<br>+  is3d = ( TYPE_HASZ(geom1->type) );
<br>+     SRID = pglwgeom_getSRID(geom1);<br>+      initGEOS(lwnotice, lwnotice);<br>+<br>+ #ifdef PGIS_DEBUG<br>+        elog(NOTICE,"simplify starting");<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+   profstart(PROF_P2G1);<br>
+ #endif<br>+     g1 = POSTGIS2GEOS(geom1);<br>+ #ifdef PROFILE<br>+  profstop(PROF_P2G1);<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE," constructed geometry - calling geos");<br>+       elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));
<br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+   profstart(PROF_GRUN);<br>+ #endif<br>+      g2 = GEOSSimplify(g1,tolerance);<br>+ #ifdef PROFILE<br>+   profstop(PROF_GRUN);
<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+   elog(NOTICE,"simplify finished");<br>+ #endif<br>+<br>+     if (g2 == NULL)<br>+      {<br>+            elog(ERROR,"GEOS Simplify() threw an error!");<br>+             GEOSGeom_destroy(g1);
<br>+             PG_RETURN_NULL(); /* never get here */<br>+       }<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+  GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE<br>+     profstart(PROF_G2P);
<br>+ #endif<br>+   result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br>+   profstop(PROF_G2P);<br>+ #endif<br>+<br>+     if (result == NULL)<br>+  {<br>+            GEOSGeom_destroy(g1);<br>+                GEOSGeom_destroy(g2);<br>+                elog(ERROR,"GEOS Simplify() threw an error (result postgis
<br>geometry formation)!");<br>+               PG_RETURN_NULL(); /* never get here */<br>+       }<br>+    GEOSGeom_destroy(g1);<br>+        GEOSGeom_destroy(g2);<br>+<br>+     /* compressType(result); */<br>+<br>+ #ifdef PROFILE<br>+     profstop(PROF_QRUN);
<br>+     profreport("geos",geom1, result);<br>+ #endif<br>+        PG_FREE_IF_COPY(geom1, 0);<br>+<br>+        PG_RETURN_POINTER(result);<br>+ }<br>+<br>+ PG_FUNCTION_INFO_V1(simplifytopologypreserve_geos);<br>+ Datum simplifytopologypreserve_geos(PG_FUNCTION_ARGS)
<br>+ {<br>+        PG_LWGEOM *geom1;<br>+    double tolerance;<br>+    GEOSGeom g1,g2;<br>+      PG_LWGEOM *result;<br>+   int is3d;<br>+    int SRID;<br>+<br>+ #ifdef PROFILE<br>+       profstart(PROF_QRUN);<br>+ #endif<br>+<br>+   geom1 = (PG_LWGEOM *)  PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
<br>+     tolerance = PG_GETARG_FLOAT8(1);<br>+<br>+  is3d = ( TYPE_HASZ(geom1->type) );<br>+     SRID = pglwgeom_getSRID(geom1);<br>+         initGEOS(lwnotice, lwnotice);<br>+<br>+ #ifdef PGIS_DEBUG<br>+        elog(NOTICE,"simplify starting");
<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+      profstart(PROF_P2G1);<br>+ #endif<br>+      g1 = POSTGIS2GEOS(geom1);<br>+ #ifdef PROFILE<br>+  profstop(PROF_P2G1);<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE," constructed geometry - calling geos");
<br>+     elog(NOTICE," g1 = %s", GEOSGeomToWKT(g1));<br>+ /*elog(NOTICE,"g1 is valid = %i",GEOSisvalid(g1)); */<br>+ #endif<br>+<br>+ #ifdef PROFILE<br>+      profstart(PROF_GRUN);<br>+ #endif<br>+      g2 = GEOSTopologyPreserveSimplify(g1,tolerance);
<br>+ #ifdef PROFILE<br>+   profstop(PROF_GRUN);<br>+ #endif<br>+<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE,"simplify finished");<br>+ #endif<br>+<br>+     if (g2 == NULL)<br>+      {<br>+            elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error!");
<br>+             GEOSGeom_destroy(g1);<br>+                PG_RETURN_NULL(); /* never get here */<br>+       }<br>+ #ifdef PGIS_DEBUG<br>+       elog(NOTICE,"result: %s", GEOSGeomToWKT(g2) ) ;<br>+ #endif<br>+  GEOSSetSRID(g2, SRID);<br>+ #ifdef PROFILE
<br>+     profstart(PROF_G2P);<br>+ #endif<br>+       result = GEOS2POSTGIS(g2, is3d);<br>+ #ifdef PROFILE<br>+   profstop(PROF_G2P);<br>+ #endif<br>+<br>+     if (result == NULL)<br>+  {<br>+            GEOSGeom_destroy(g1);<br>+                GEOSGeom_destroy(g2);
<br>+             elog(ERROR,"GEOS TopologyPreserveSimplify() threw an error (result<br>postgis geometry formation)!");<br>+                PG_RETURN_NULL(); /* never get here */<br>+       }<br>+    GEOSGeom_destroy(g1);<br>+        GEOSGeom_destroy(g2);
<br>+<br>+  /* compressType(result); */<br>+<br>+ #ifdef PROFILE<br>+     profstop(PROF_QRUN);<br>+         profreport("geos",geom1, result);<br>+ #endif<br>+        PG_FREE_IF_COPY(geom1, 0);<br>+<br>+        PG_RETURN_POINTER(result);
<br>+ }<br>+<br>diff -r -c /postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp<br>/new/postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp<br>*** /postgis-1.1.6/lwgeom/lwgeom_geos_wrapper.cpp        Fri Dec 16 00:56:15 2005<br>--- /new/postgis-
1.1.6/lwgeom/lwgeom_geos_wrapper.cpp    Fri Nov 17 07:54:29 2006<br>***************<br>*** 10,15 ****<br>--- 10,21 ----<br>  #include "postgis_geos_version.h"<br>   #include "geos/geom.h"<br>   #include "geos/util.h"
<br>+<br>+<br>+ //#include <geos/simplify/DouglasPeuckerSimplifier.h><br>+ //#include <geos/simplify/TopologyPreservingSimplifier.h><br>+<br>+<br>  #if GEOS_FIRST_INTERFACE <= 3 && GEOS_LAST_INTERFACE >= 3
<br>   #include "geos/opValid.h"<br>   #include "geos/opPolygonize.h"<br>***************<br>*** 140,145 ****<br>--- 146,155 ----<br>   extern "C" Geometry *GEOSSymDifference(Geometry *g1,Geometry *g2);
<br>   extern "C" Geometry *GEOSUnion(Geometry *g1,Geometry *g2);<br><br>+<br>+ extern "C" Geometry *GEOSSimplify(Geometry *g1, double tolerance);<br>+ extern "C" Geometry *GEOSTopologyPreserveSimplify(Geometry *g1,
<br>double tolerance);<br>+<br><br>   extern "C" POINT3D  *GEOSGetCoordinate(Geometry *g1);<br>   extern "C" POINT3D  *GEOSGetCoordinates_Polygon(Polygon *g1);<br>***************<br>*** 1857,1859 ****<br>
--- 1867,1917 ----<br>            return NULL;<br>          }<br>  }<br>+<br>+<br>+ Geometry *<br>+ GEOSSimplify(const Geometry *g1, double tolerance)<br>+ {<br>+        using namespace geos::simplify;<br>+<br>+   try<br>+  {<br>+            Geometry::AutoPtr g(DouglasPeuckerSimplifier::simplify(g1, tolerance));
<br>+             return g.release();<br>+  }<br>+    catch (const std::exception &e)<br>+  {<br>+            ERROR_MESSAGE("%s", e.what());<br>+             return NULL;<br>+         }<br>+<br>+         catch (...)<br>+  {<br>+            ERROR_MESSAGE("Unknown exception thrown");
<br>+             return NULL;<br>+         }<br>+ }<br>+<br>+ Geometry *<br>+ GEOSTopologyPreserveSimplify(const Geometry *g1, double tolerance)<br>+ {<br>+   using namespace geos::simplify;<br>+<br>+   try<br>+  {<br>+            Geometry::AutoPtr g(TopologyPreservingSimplifier::simplify(g1, tolerance));
<br>+             return g.release();<br>+  }<br>+    catch (const std::exception &e)<br>+  {<br>+            ERROR_MESSAGE("%s", e.what());<br>+             return NULL;<br>+         }<br>+<br>+         catch (...)<br>+  {<br>+            ERROR_MESSAGE("Unknown exception thrown");
<br>+             return NULL;<br>+         }<br>+ }<br>+<br>diff -r -c /postgis-1.1.6/lwgeom/lwpostgis.sql.in<br>/new/postgis-1.1.6/lwgeom/lwpostgis.sql.in<br>*** /postgis-1.1.6/lwgeom/lwpostgis.sql.in    Fri Jul  7 06:56:52 2006<br>--- /new/postgis-
1.1.6/lwgeom/lwpostgis.sql.in   Fri Nov 17 08:55:30 2006<br>***************<br>*** 2876,2881 ****<br>--- 2876,2893 ----<br>   -- GEOS<br>  ---------------------------------------------------------------<br><br>+<br>+ CREATEFUNCTION simplify_geos(geometry,float8)
<br>+    RETURNS geometry<br>+    AS '@MODULE_FILENAME@','simplify_geos'<br>+    LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable)<br>+<br>+ CREATEFUNCTION simplifytopologypreserve_geos(geometry,float8)<br>+    RETURNS geometry
<br>+    AS '@MODULE_FILENAME@','simplifytopologypreserve_geos'<br>+    LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable)<br>+<br>+<br>   CREATEFUNCTION intersection(geometry,geometry)<br>     RETURNS geometry
<br>     AS '@MODULE_FILENAME@','intersection'</pre>