Re: [postgis-users] Speeding up point-in-polygon search using ST_SimplifyPreserveTopology?

[Evan Martin]

Thanks, Paul, but do you have any suggestions on how I can cut up the 
polygons properly? As I posted, ST_Intersection doesn't work properly on 
the ones spanning the dateline and I don't know how else to do it. Maybe 
there's some obvious trick I'm missing. (I haven't looked at the code, 
but it seems a bit strange to me that ST_Intersects works, but not 
ST_Intersection - so PostGIS can somehow figure out that they intersect, 
but can't figure out where?)


By the development work do you mean ticket 1796 or something else? 
Either way, any estimate on when it might get done? :)


Regards,

Evan

On 30/05/2012 11:35 PM, Paul Ramsey wrote:

Evan,
Beyond the cutting-up-of-large-things as Stephen suggests, there is
the extra development work to optimize the edge searching routines
with an internal index on the objects (and caching that index once
built). Something for 2.1 :)
P.

On Mon, May 28, 2012 at 10:08 AM, Evan Martin
postgre...@realityexists.net  wrote:

Hi all,

I have a table of ~350 multi-polygons and ~185,000 points and I need to find
which points are inside which multi-polygons. Some polygons are quite large
and span the dateline, so I'm using geography ST_DWithin for this (with a
tolerance of 100m). My initial query looks like this (simplified):

SELECT ...
FROM points, polygons
WHERE ST_DWithin(point, real_area, 100)

This works, but it takes about 90 minutes. I'm trying to make it faster by
using ST_SimplifyPreserveTopology. That worked very nicely for my adjacent
polygons problem
[http://postgis.refractions.net/pipermail/postgis-users/2012-January/031992.html],
because all polygons were modified in the same way, but this is trickier.
Since I'm modifying the polygon, but not the point, the results are
different. So I thought maybe I could do this in two phases: if the point is
well within or well outside the polygon then take the result of the
simplified check as correct, but if it's close to the border then check it
properly, ie.

SELECT ...
FROM points, polygons
WHERE ST_DWithin(point, simple_area, 2)
AND (ST_Distance(point, simple_border)  2 OR ST_DWithin(point,
real_area, 100))

simple_area = ST_SimplifyPreserveTopology(real_area::geometry,
0.01)::geography and simple_border = the exterior ring of simple_area.

This takes about 18 minutes (a 5x improvement) and gives very similar
results, but not quite the same. It falls down on polygons that have
rhumblines along parallels, because they get turned into great circle lines.
Eg. the original polyon may have a rhumbline approximated as (24 10,25 10,26
10,27 10), ST_SimplifyPreserveTopology does its job and simplifies it to (24
10,27 10) and then ST_DWithin on geography treats it as a great circle line,
giving an incorrect result. I tried inserting extra points to unsimplify
the rhumblines, but that itself is very slow and also quite a hack, because
I can't really be sure which lines were supposed be rhumblines when looking
at the simplified polygon. I feel like I'm so close and this is a very silly
little problem - but it is a problem.

Could anyone suggest a way to work around this? Or a different approach to
the whole thing?

Thanks,

Evan
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users


___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

Re: [postgis-users] Speeding up point-in-polygon search using ST_SimplifyPreserveTopology?

[Evan Martin]

Thanks, Stephen, but I haven't been able to figure out how to get this 
to work on geography. ST_Intersection on geography doesn't work properly 
across the dateline. For example, I tried manually tiling a simple 
polygon that goes across the dateline, like this:


WITH data AS
(
SELECT ST_GeomFromEWKT('SRID=4326;POLYGON((163 50,-176 55,-176 
60,163 60,163 50))')::geography AS simple

),
tiles AS
(
SELECT ST_SetSRID('BOX(145 50,150 60)'::box2d, 4326) AS tile
UNION ALL SELECT ST_SetSRID('BOX(150 50,155 60)'::box2d, 4326)
UNION ALL SELECT ST_SetSRID('BOX(155 50,160 60)'::box2d, 4326)
UNION ALL SELECT ST_SetSRID('BOX(160 50,165 60)'::box2d, 4326)
UNION ALL SELECT ST_SetSRID('BOX(165 50,170 60)'::box2d, 4326)
UNION ALL SELECT ST_SetSRID('BOX(170 50,175 60)'::box2d, 4326)
UNION ALL SELECT ST_SetSRID('BOX(175 50,180 60)'::box2d, 4326)
UNION ALL SELECT ST_SetSRID('BOX(-180 50,-175 60)'::box2d, 4326)
UNION ALL SELECT ST_SetSRID('BOX(-175 50,-170 60)'::box2d, 4326)
)
SELECT ST_AsText(ST_Intersection(simple, tile::geography)) AS tiled
FROM data
JOIN tiles ON ST_Intersects(simple, tile::geography)

This returns:

POLYGON((163 49.99887,160 50.0467770263591,160 
59.99987,163 59.99987,163 49.99887))

GEOMETRYCOLLECTION EMPTY
GEOMETRYCOLLECTION EMPTY
GEOMETRYCOLLECTION EMPTY
POLYGON((-175 54.9860851802267,-176 54.99957,-176 
59.99987,-175 59.99987,-175 54.9860851802267))


I think ST_Intersects works OK on geography, but not ST_Intersection.

Regards,

Evan


On 29/05/2012 2:21 AM, Stephen Woodbridge wrote:

On 5/28/2012 10:08 AM, Evan Martin wrote:

Hi all,

I have a table of ~350 multi-polygons and ~185,000 points and I need to
find which points are inside which multi-polygons. Some polygons are
quite large and span the dateline, so I'm using geography ST_DWithin for
this (with a tolerance of 100m). My initial query looks like this
(simplified):

SELECT ...
FROM points, polygons
WHERE ST_DWithin(point, real_area, 100)

This works, but it takes about 90 minutes. I'm trying to make it faster
by using ST_SimplifyPreserveTopology. That worked very nicely for my
adjacent polygons problem
[http://postgis.refractions.net/pipermail/postgis-users/2012-January/031992.html], 


because all polygons were modified in the same way, but this is
trickier. Since I'm modifying the polygon, but not the point, the
results are different. So I thought maybe I could do this in two phases:
if the point is well within or well outside the polygon then take the
result of the simplified check as correct, but if it's close to the
border then check it properly, ie.

SELECT ...
FROM points, polygons
WHERE ST_DWithin(point, simple_area, 2)
AND (ST_Distance(point, simple_border)  2 OR ST_DWithin(point,
real_area, 100))

simple_area = ST_SimplifyPreserveTopology(real_area::geometry,
0.01)::geography and simple_border = the exterior ring of simple_area.

This takes about 18 minutes (a 5x improvement) and gives very similar
results, but not quite the same. It falls down on polygons that have
rhumblines along parallels, because they get turned into great circle
lines. Eg. the original polyon may have a rhumbline approximated as (24
10,25 10,26 10,27 10), ST_SimplifyPreserveTopology does its job and
simplifies it to (24 10,27 10) and then ST_DWithin on geography treats
it as a great circle line, giving an incorrect result. I tried inserting
extra points to unsimplify the rhumblines, but that itself is very
slow and also quite a hack, because I can't really be sure which lines
were supposed be rhumblines when looking at the simplified polygon. I
feel like I'm so close and this is a very silly little problem - but it
is a problem.

Could anyone suggest a way to work around this? Or a different approach
to the whole thing?


I think the alternative approach that has been discussed on the list 
in the past, and should be in the archives, is to cut the multiple 
polygons into tiles with all the attributes of the original 
multipolygon and then to match the points to the tiles..


This works much faster for two basic reasons:

1. the number of points in the each tile is much less than the 
original because it only contains a fraction of the originals complex 
boundary but maintains all the detail.
2. since each tile is spatially smaller, you get better (faster) 
interactions between the tile index and the points.


-Steve
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users


___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[postgis-users] Speeding up point-in-polygon search using ST_SimplifyPreserveTopology?

[Evan Martin]

Hi all,

I have a table of ~350 multi-polygons and ~185,000 points and I need to 
find which points are inside which multi-polygons. Some polygons are 
quite large and span the dateline, so I'm using geography ST_DWithin for 
this (with a tolerance of 100m). My initial query looks like this 
(simplified):


SELECT ...
FROM points, polygons
WHERE ST_DWithin(point, real_area, 100)

This works, but it takes about 90 minutes. I'm trying to make it faster 
by using ST_SimplifyPreserveTopology. That worked very nicely for my 
adjacent polygons problem 
[http://postgis.refractions.net/pipermail/postgis-users/2012-January/031992.html], 
because all polygons were modified in the same way, but this is 
trickier. Since I'm modifying the polygon, but not the point, the 
results are different. So I thought maybe I could do this in two phases: 
if the point is well within or well outside the polygon then take the 
result of the simplified check as correct, but if it's close to the 
border then check it properly, ie.


SELECT ...
FROM points, polygons
WHERE ST_DWithin(point, simple_area, 2)
AND (ST_Distance(point, simple_border)  2 OR ST_DWithin(point, 
real_area, 100))


simple_area = ST_SimplifyPreserveTopology(real_area::geometry, 
0.01)::geography and simple_border = the exterior ring of simple_area.


This takes about 18 minutes (a 5x improvement) and gives very similar 
results, but not quite the same. It falls down on polygons that have 
rhumblines along parallels, because they get turned into great circle 
lines. Eg. the original polyon may have a rhumbline approximated as (24 
10,25 10,26 10,27 10), ST_SimplifyPreserveTopology does its job and 
simplifies it to (24 10,27 10) and then ST_DWithin on geography treats 
it as a great circle line, giving an incorrect result. I tried inserting 
extra points to unsimplify the rhumblines, but that itself is very 
slow and also quite a hack, because I can't really be sure which lines 
were supposed be rhumblines when looking at the simplified polygon. I 
feel like I'm so close and this is a very silly little problem - but it 
is a problem.


Could anyone suggest a way to work around this? Or a different approach 
to the whole thing?


Thanks,

Evan
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

Re: [postgis-users] Transform from NAD27 (SRID 4267) to WGS84 (SRID 4326)

[Evan Martin]

70 west, 10 north is in venezuela which is outside the bounds of
available NAD27 to NAD83 conversion files.  Try (-120 40).

Thanks, Frank. Yes, (-120 40) works. Just to confirm, will PostGIS 
always return the input coordinates now when they're outside the bounds 
of the input datum, rather than throwing an error? Is that the expected 
behaviour?


Regards,

Evan
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[postgis-users] Transform from NAD27 (SRID 4267) to WGS84 (SRID 4326)

[Evan Martin]

When trying to transform from NAD27 (SRID 4267) to WGS84 (SRID 4326) 
some points used to fail as of 2.0.0 beta 3 (bug 318). I've just tried 
2.0.0 RC 2 and they now return the input value unchanged, eg.


SELECT ST_AsText(ST_Transform(ST_GeomFromEWKT('SRID=4267;POINT(-70 
10)'), 4326))


returns POINT(-70 10). That's definitely better than an error, but is 
that result correct? I thought there was some difference between NAD27 
and WGS84. CoordTrans (http://franson.com/coordtrans/) returns 
(-69.9995694 10.001074) for the same transformation.


Thanks for any help,

Evan
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[postgis-users] Find polygon border crossings over dateline

[Evan Martin]

I have a table of polygons covering most of the world. I posted a 
question before about finding which of them are adjacent: 
http://postgis.refractions.net/pipermail/postgis-users/2012-January/031992.html 
Now I also have a set of great-circle lines (LINESTRINGs with only 2 
points) and I need to figure out which of them cross the border between 
two given polygons (in the right direction).


As before, there is often no precise border between the polygons - 
sometimes they overlap a little and sometimes there is a slight gap. So 
it's not just a matter of calling ST_Intersects(border_line, 
crossing_line). I tried 2 ways of doing this:


1) Buffer each polygon a bit to fill in the gaps, get the intersections 
and check if the line crosses that, ie.
ST_Intersects(ST_Intersection(ST_Buffer(polygon1, 100), 
ST_Buffer(polygon2, 100)), crossing_line)


2) Find ST_Intersection(polygon, crossing_line) for each polygon, which 
is a LINESTRING with 2 points. Where the end point of one of these 
LINESTRINGs is the start point of another that means crossing_line went 
from polygon1 to polygon2.


Both of these approaches seem to basically work - except when the shapes 
span the dateline. Doing the calculations on geometry doesn't return the 
right results, because the coordinates differ between shapes (eg. 190 
vs. -170). Doing them on geography fails with this error:


ERROR:  Error performing intersection: TopologyException: found 
non-noded intersection between LINESTRING ...


I'm guessing that's because ST_Buffer and ST_Intersection on geography 
really work on geometry internally anyway.


Any ideas on how I can do this? Thanks in advance for any help!

Evan
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[postgis-users] Normalise geometry for conversion to geography

[Evan Martin]

Is there any built-in function or other simple way to normalise a 
geometry to the (-180-180, -90-90) range required by geography? I've 
written my own function for this, but it's quite slow, because I need to 
handle MultiPolygons so if even a single point needs to be changed it 
has to re-create the whole thing point-by-point.


I found this discussion about doing it automatically when casting to 
geography: 
http://postgis.refractions.net/pipermail/postgis-devel/2011-June/013805.html 
and I can see why that might be controversial, but what about some way 
for the user to explicitly request normalization?


Thanks,

Evan
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[postgis-users] A faster way to find adjacent polygons

[Evan Martin]

I have a table of polygons covering most of the world and for a given 
polygon I need to find all polygons adjacent to it. In theory, two 
adjacent polygons should have a part of their border in common, but in 
practice the data is not that precise - sometimes there's a slight 
overlap and sometimes a slight gap between adjacent polygons. Since 
these polygons can cross the dateline I have them stored as geography. 
It looks like ST_DWithin(geog, geog, distance) does exactly what I want, eg.


SELECT *
FROM polygons p1
JOIN polygons p2 ON p1.id  p2.id AND ST_DWithin(p1.geog, p2.geog, 100, 
false)


The only problem is that this is slow! There are about 300 polygons in 
total and this query takes about 75 minutes. I have a GIST index defined 
on the geography column. Just doing  is fast - it seems to be 
_ST_DWithin that's slow. There's one pair of polygons I found (with 300 
points and 1300 points) that share a part of the border (ST_Intersects = 
true) and _ST_DWithin takes 1.2 or 3.9 seconds just for that one pair, 
depending on the order of the arguments.


Is there some other way I can do this? I couldn't figure out any way to 
get the correct result with geometry around the dateline and ST_DWithin 
looks like the appropriate geography function.


Can ST_DWithin itself be made faster? I had a look at the source and it 
seems to be calculating the distance between each combination of edges - 
so O(n^2). Isn't there a more efficient algorithm for this out there? I 
know there's Rotating Calipers for a plane - is there something like 
that which works on a sphere? Or is it possible to consider only the 
edges of each polygon that are within the intersection of their bounding 
boxes? I'm a newbie to this, so it might be a silly idea, but surely 
there has to be something faster than the brute force approach?


Any help would be appreciated!

Evan

___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[postgis-users] ST_SnapToGrid returns a value out of range

[Evan Martin]

I've found some odd behaviour in ST_SnapToGrid, which looks like a bug. 
When I run


SELECT ST_SnapToGrid(ST_Transform(ST_GeomFromText('POINT(180 50)', 
4269), 4326), 0.1)::geography


an error is returned:

ERROR:  Coordinate values are out of range [-180 -90, 180 90] for 
GEOGRAPHY type


Calling ST_AsText() on the same geometry returns POINT(180 50), as 
expected. Interestingly, if I either decrease or even increase the 
precision the error does not occur. Eg. this works fine:


SELECT ST_SnapToGrid(ST_Transform(ST_GeomFromText('POINT(180 50)', 
4269), 4326), 0.01)::geography


It's only the value 0.0001 that seems to cause a problem. Is this a bug 
or am I missing something?


Regards,

Evan
___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users

[postgis-users] Geography and CIRCULARSTRING / CURVEPOLYGON

[Evan Martin]

I'm using PostGIS 2.0.0 and would like to store a GML Surface in a way 
that allows calculating intersections between that surface and various 
points and lines. In my case, the Surface will have one PolygonPatch 
with no interior rings and an exterior ring which may be a 
LineStringSegment, a Geodesic, an Arc or a Circle. The distances can be 
in the order of tens to hundreds of miles.


CIRCULARSTRING and CURVEPOLYGON would seem like a good way of storing an 
Arc, except that they're not supported by the geography type, only 
geometry. Are there any plans to add support for them any time soon? 
What would be the best approach to this in the meanwhile?


Regards,

Evan Martin

___
postgis-users mailing list
postgis-users@postgis.refractions.net
http://postgis.refractions.net/mailman/listinfo/postgis-users