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
_______________________________________________
postgis-users mailing list
[email protected]
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
