Hi
Thanks. (Sorry for the late replay, we where upgrading out main database in DMZ
from Postgres 9.6 to Postgres 11.3 so I had to work with that for a couple of
days.)
Now I have done some testing comparing srid 3035 and said 4258. I stared out
with a dataset representing all the municipalities in Norway. This is a
dataset with no overlaps.
So I created a new dataset using srid 3035
create table test_srid_3035.data_as_srid_3035 as (select sl_sdeid,
ST_transform(geo,3035)::geometry(Polygon,3035) as geom from
sk_grl.n5_forvaltning_flate);
This new data did not have any gaps either, so that was very good.
Then I transformed it back 4258
create table test_srid_3035.data_as_srid_3035_to_4258 as (select sl_sdeid,
ST_transform(geom,4258)::geometry(Polygon,4258) as geom from
test_srid_3035.data_as_srid_3035);
and it was still perfect.
Then I made two new data with points based on the municipality dataset.
create table test_srid_3035.test_4258_point as (select sl_sdeid,
ST_Centroid(geo)::geometry(Point,4258) as geom from
sk_grl.n5_forvaltning_flate);
create table test_srid_3035.test_3035_point as (select sl_sdeid,
ST_transform(ST_Centroid(geo),3035)::geometry(Point,3035) as geom from
sk_grl.n5_forvaltning_flate);
And I did your distance test like you suggested on the 4258 table . The total
numbers of distanses testet on is 182756.
select min(distance) as min_4258_spheroid_true, max (distance) as
max_4258_spheroid_true, avg(distance) as avg_4258_spheroid_true, sum(distance)
as sum_4258_spheroid_true from
(
select ST_Distance(t1.geom::geography,t2.geom::geography) as distance
from
test_srid_3035.test_4258_point t1,
test_srid_3035.test_4258_point t2
where t1.sl_sdeid != t2.sl_sdeid
order by t1.sl_sdeid, t2.sl_sdeid
) as r;
min_4258_spheroid_true | max_4258_spheroid_true | avg_4258_spheroid_true |
sum_4258_spheroid_true
------------------------+------------------------+------------------------+------------------------
731.66154569 | 1787695.9610192 | 507972.286659922 |
92834983220.8208
(1 row)
Time: 362.910 ms
The I tested 3035 table like this.
select min(distance) as min_3035, max (distance) as max_3035, avg(distance) as
avg_3035, sum(distance) as sum_3035 from
(
select ST_Distance(t1.geom,t2.geom) as distance
from
test_srid_3035.test_3035_point t1,
test_srid_3035.test_3035_point t2
where t1.sl_sdeid != t2.sl_sdeid
order by t1.sl_sdeid, t2.sl_sdeid
) as r;
min_3035 | max_3035 | avg_3035 | sum_3035
------------------+------------------+------------------+-----------------
735.136469548696 | 1777644.86377117 | 506075.903269244 | 92488407777.874
(1 row)
Time: 145.513
We see differences here as probably expected .
I did one more distance test where transferred the points back to 4258 from
3035.
select min(distance) as min_3035_backto_4258, max (distance) as
max_3035_backto_4258, avg(distance) as avg_3035_backto_4258, sum(distance) as
sum_3035_backto_4258 from
(
select
ST_Distance(ST_transform(t1.geom,4258)::geography,ST_transform(t2.geom,4258)::geography)
as distance
from
test_srid_3035.test_3035_point t1,
test_srid_3035.test_3035_point t2
where t1.sl_sdeid != t2.sl_sdeid
order by t1.sl_sdeid, t2.sl_sdeid
) as r;
min_3035_backto_4258 | max_3035_backto_4258 | avg_3035_backto_4258 |
sum_3035_backto_4258
----------------------+----------------------+----------------------+----------------------
731.66154572 | 1787695.96111941 | 507972.286689737 |
92834983226.2696
(1 row)
Time: 498.691 ms
Then we get some very,verry minor differences.
So this 3035 looks quite ok , if your main focus is not to work on distances
between objects. If you need exact correct distance you have to transform it to
geography .
So is there any good reasons not to use 3035 if your main focus is the keep
track of polygon borders and the areas of polygon ?
Thanks
Lars
________________________________
From: Paul Ramsey <[email protected]>
Sent: Tuesday, June 11, 2019 5:36 PM
To: Lars Aksel Opsahl
Cc: PostGIS Users Discussion; Sandro Santilli
Subject: Re: [postgis-users] Tolerance/SnapTo in Postgis Topology for meter and
degrees.
Not ideal, a whole-of-Norway projection would have even more error
minimization, but not a bad stopgap, certainly better than using UTM
projections way outside their areas of validity. Check a few measurements, and
do things like comparing an ST_Distance(A, B) with ST_Distance(A::geography,
B::geography) to see how much deviation from truth the projection is adding.
P.
On Jun 11, 2019, at 4:34 AM, Lars Aksel Opsahl
<[email protected]<mailto:[email protected]>> wrote:
Hi
I had a discussion with Knut Bjørkelo a co worker of me today.
He was wondering if we could use https://epsg.io/3035, then all the area values
should be correct for all valid UTM Zones in Noway.
When we need to show a map we transform in to the correct to UTM zone to get a
visual correct map also.
If I test a single point the transformation seems be accurate up to 4. decimal
using 3035 even for a extreme point.
select
ST_AStexT(ST_Transform(ST_Transform(ST_Transform(ST_Transform(ST_setSrid(ST_MakePoint(1108142.0,7788000.0),25835),3035),25832),3035),25835));
st_astext
------------------------------------------
POINT(1108142.00005503 7787999.99979103)
(1 row)
Lars
________________________________
From: postgis-users
<[email protected]<mailto:[email protected]>>
on behalf of Lars Aksel Opsahl
<[email protected]<mailto:[email protected]>>
Sent: Monday, June 10, 2019 10:50 PM
To: PostGIS Users Discussion
Cc: Sandro Santilli
Subject: Re: [postgis-users] Tolerance/SnapTo in Postgis Topology for meter and
degrees.
Hi
You have this C-code
4870 _lwt_minTolerance( LWGEOM *g )
4871 {
4872 const GBOX* gbox;
4873 double max;
4874 double ret;
4875
4876 gbox = lwgeom_get_bbox(g);
4877 if ( ! gbox ) return 0; /* empty */
4878 max = FP_ABS(gbox->xmin);
4879 if ( max < FP_ABS(gbox->xmax) ) max = FP_ABS(gbox->xmax);
4880 if ( max < FP_ABS(gbox->ymin) ) max = FP_ABS(gbox->ymin);
4881 if ( max < FP_ABS(gbox->ymax) ) max = FP_ABS(gbox->ymax);
4882
4883 ret = 3.6 * pow(10, - ( 15 - log10(max?max:1.0) ) );
4884
4885 return ret;
4886 }
In postgis branch svn-2.5 we have this SQL code.
CREATE OR REPLACE FUNCTION topology._st_mintolerance(ageom Geometry)
RETURNS float8
AS $$
SELECT 3.6 * power(10, - ( 15 - log(coalesce(
nullif(
greatest(abs(ST_xmin($1)), abs(ST_ymin($1)),
abs(ST_xmax($1)), abs(ST_ymax($1))),
0),
1)) ));
$$ LANGUAGE 'sql' IMMUTABLE STRICT;
I am not that into C-code but the SQL code looks like the C-code above as you
suggested in your mail.
So I did a small test
https://github.com/NibioOpenSource/pgtopo_update_sql/blob/develop/src/test/sql/snapto/snapto_code_example_degrees_st_min_tolerance.sql
using topology._st_mintolerance and yes I was able to get snap to every second
line both vertically and horizontally to work and both and north and south in
Norway. I just got this to work by playing around with
topology._st_mintolerance together different factor values.
But I have no idea about how I can use this min tolerance code to make a
generic code that behaves like 10 meter tolerance in a planar projection.
And I am also little unsure about what projection that is correct to use for
datasets that covers all of Norway now.
Before we did not really have any choice because the projection was quite bad
if moved outside the bonds of the projection limits. We pretty much had to use
degrees to get the a OK quality if we needed to store the data in one single
dataset.
So whats the "correct projection" now to use for original dataset that covers
all of Norway and which are updated regularly ?
Thanks a lot.
Lars
________________________________
From: [email protected]<mailto:[email protected]> <[email protected]<mailto:[email protected]>> on
behalf of Sandro Santilli <[email protected]<mailto:[email protected]>>
Sent: Saturday, June 8, 2019 2:06 PM
To: PostGIS Users Discussion
Cc: Lars Aksel Opsahl
Subject: Re: [postgis-users] Tolerance/SnapTo in Postgis Topology for meter and
degrees.
On Fri, Jun 07, 2019 at 10:02:23AM +0000, Lars Aksel Opsahl wrote:
> So the problem is how to use tolerances so we get a behavior equal to
> the test using meter.
>
>
> We can we define the layer in Postgis Topology with quite big value
> because this is just max value as it seems. So we can adjust the tolerance
> parameter as we add lines but the problem is that we need to adjust this
> parameter depending on where we are and what orientation the line has. For
> vertical lines we need a bigger tolerance than for horizontal lines in
> Norway. This makes it quite complicated to handle adding new lines.
I'm not sure if it'd help but PostGIS Topology has an internal
function (not to be relied upon, but could be copied to your own
function) to determine "min tolerance" based on absolute coordinate
values. That function is meant to deal with non-uniform floating-point
resolution. What you're after is a function to deal with non-uniform
tolerances. See
https://trac.osgeo.org/postgis/browser/trunk/liblwgeom/lwgeom_topo.c?rev=14251#L4866
It used to be done in SQL with previous versions.
Could that help ?
--strk;
_______________________________________________
postgis-users mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/postgis-users
