[postgis-users] ST_Simplify and projections

[Paolo Crosato]

Il 08/01/2013 13:42, Sandro Santilli ha scritto:
> Hi Paolo,
>
> On Tue, Jan 08, 2013 at 01:02:15PM +0100, Paolo Crosato wrote:
>> Hello,
>>
>> I'm trying to simplify a set of administrative boundaries that spans
>> over the whole northern emishpere, boundaries are:
>>
>> xMin,yMin -180,27.6378 : xMax,yMax 180,83.6274, CRS is 4326 WGS84.
>> The whole dataset comprises 100k features and 54 million vertices,
>> so it's quite big.
>>
>>
>> Since I have to preserve the topology, I'm using topologies to do
>> the simplification, as in this article: http://strk.keybit.net/blog/2012/04/13/simplifying-a-map-layer-using-postgis-topology/.
>> The simplification process uses ST_simplify, and I know that using
>> this function with non planar projections is not recommended. I
>> still get good results with 4326, even at medium scales.
>> However at large scales, like 1:17893297 or so (zoom level 5 to 1 on
>> a tiled map system), there are visible artifacts. I've tried to
>> simplify using lambert projection
>> (EPSG:2154) on a single country, and the results are slightly better.
>>
>>
>> I'm not an expert of geographical projections, so I'd like to know
>> if there is some sort of planar projection supported by postgis,
>> that spans over the whole northern emisphere.
>> Visual distortion is not an issue, my idea is to project the whole
>> set to this new CRS, simplify and reproject back the results to
>> 4326. Simplifying single countries of the dataset using different
>> projections is not feasible, since I have to simplify the whole set
>> at once to preserve the topology.
> Why not simplify single _edges_ of the topology using different
> projections ? I mean, for each edge:
>   - project to best projection for it
>   - simplify
>   - project back
> Doing the above in the loop performed by my wrapper function (in the
> blog post) should guard against topology errors.
>
> SRID of "best projection" may currently be done using the internal
> and undocumented _ST_BestSRID(geography) function.
> It may actually be interesting to expose an
> ST_Simplify(geography, <meters>) to make this even simpler.
>
Hi,

thank you very much for the suggestion, didn't thought about it. I 
wrapped all the code in a SimplifyUTM function, and I call it from your 
function instead of simplify.
I preferred the usual utmzone() function, which is found in the postgis 
book, instead of the _ST_BestSRID function. Tests on the prototype were 
very good, I hope the utm
transformation doesn't give any problem in Siberia or Kazakhstan, where 
administrative boundaries have very long edges. But I need to test the 
function against
the whole european set in order to have any answer, and it's still 
building the topology :)
Anyway here is the function, maybe someone else could use it:

CREATE OR REPLACE FUNCTION simplifyutm(geomin geometry, meters double 
precision)
   RETURNS geometry AS
$BODY$
DECLARE
     bestUTM integer;
     geomUTM geometry;
     geomSimpl geometry;
     geomBack geometry;
BEGIN
     bestUTM:= utmzone(st_centroid(geomin));
     geomUTM:= ST_Transform(geomin,bestUTM);
     geomSimpl:= ST_Simplify(geomUTM,meters);
     geomBack:= ST_Transform(geomSimpl,ST_SRID(geomin));
     geomBack:= ST_SnapToGrid(geomBack,1e-5);
     geomBack:= ST_SetPoint(geomBack, 0, ST_StartPoint(geomin));
     geomBack:= ST_SetPoint(geomBack, ST_NumPoints(geomBack)-1, 
ST_EndPoint(geomin));
     RETURN geomBack;
END
$BODY$
   LANGUAGE plpgsql STABLE STRICT;

Caveats:
I'm a newbie at pg/SQL functions, so this could be better rewritten, at 
least to handle exceptions.
It's tailored to my needs, so it discards any decimal digit beyond the 
5th and replace start and end point with the input's ones, to preserve 
topology, so it's working on linear geometries only.

Regards,

Paolo

-- 
Paolo Crosato