<div dir="ltr"><div><div><div><div><div><div><div><div><div>Hey,<br></div><div>@Paul<br></div>Paul I'm afraid there is a true postgis problem here.<br></div>ST_Intersects for geography amount to <br>`( geom1 && geom2 AND _ST_Distance(geom1, geom2 , 0.0, false) < 0.00001 ) `<br></div>now manually checking, && gives correct result, <b>but `_ST_Distance` returns wrong results in this case</b>,<br></div><br>Yet `_st_distanceuncached(geom, multip, FALSE)` <b>returns correct answer (0).</b><br></div>I checked also by putting "0." before each coordinate to rule out a too-big-value issue.<br></div><div>I also ruled out a mixed srid problem.<br></div><div>At the end I put the complete test script, I must stop here.<br></div><div><br></div><div>@Bob<br></div>So Bob for the moment, you have some work around :<br></div> * use geometry type if you can afford<br></div> * replace all "ST_Intersects(geom1,geom2)" by exactly this :<br> ( (geom1&& geom2) AND ( _st_distanceuncached(geom1, geom2, FALSE) < 0.00001) ) <br></div><div> * replace all "ST_Intersects(geom1,geom2)" by exactly this :<br></div><div> ST_IsEmpty(ST_Intersection(geom1,geom2)::geometry)=FALSE<br><br></div><div>I don't know which one of option 2 or 3 will be slowest. It depends on your data.<br></div>Note that it will be slow any way if you have many geometry that have intersecting bounding boxes.<br><div><div><div><div><div><div><div><div><div><div><div><div><br>Cheers,<br>Rémi-C<br><br><br>----------------------<br>CREATE SCHEMA errror_distance_geography;<br>SET search_path to errror_distance_geography, public ; <br>DROP TABLE IF EXISTS multipolygon ;<br>create table multipolygon(<br> mpid SERIAL not null,<br> multip Geography(multipolygon) not null<br>);<br>insert into multipolygon <br>values (1, ST_GeographyFromText('multipolygon(((0.20 0.10, 0.21 0.10, 0.22 0.10, 0.23 0.10, 0.24 0.10, 0.25 0.10, 0.26 0.10, 0.27 0.10, 0.28 0.10, 0.29 0.10, 0.30 0.10, 0.30 0.20, 0.29 0.20, 0.28 0.20, 0.27 0.20, 0.26 0.20, 0.25 0.20, 0.24 0.20, 0.23 0.20, 0.22 0.20, 0.21 0.20, 0.20 0.20, 0.20 0.10)))'));<br>insert into multipolygon <br>values (2, ST_GeographyFromText('multipolygon(((0.25 0.15, 0.26 0.15, 0.27 0.15, 0.28 0.15, 0.29 0.15, 0.30 0.15, 0.31 0.15, 0.32 0.15, 0.33 0.15, 0.34 0.15, 0.35 0.15, 0.35 0.25, 0.34 0.25, 0.33 0.25, 0.32 0.25, 0.31 0.25, 0.30 0.25, 0.29 0.25, 0.28 0.25, 0.27 0.25, 0.26 0.25, 0.25 0.25, 0.25 0.15)))'));<br>SELECT *<br>FROM multipolygon ; <br>WITH the_geom_to_intersect AS (<br> SELECT ST_GeographyFromText('MULTIPOLYGON(((0.27 0.17, 0.28 0.17, 0.28 0.19, 0.27 0.19, 0.27 0.17)))') AS geom<br>) <br>SELECT <br> ( geom && multip AND _ST_Distance(geom, multip , 0.0, false) < 0.00001 ) AS intersects_for_geography<br> , _ST_Distance(geom,multip, 0.0, false) AS distance_for_geography<br> ,_st_distanceuncached(geom, multip, FALSE) AS uncached_distance_for_geography<br> ,( geom && multip AND _st_distanceuncached(geom, multip, FALSE) < 0.00001 ) intersects_for_geography_uncached<br> ,ST_IsEmpty(ST_Intersection(geom,multip)::geometry) is_intersection_empty_trick<br>from the_geom_to_intersect, multipolygon m1 <br>---------------------<br><br></div></div></div></div></div></div></div></div></div></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-01-09 0:44 GMT+01:00 Paul Ramsey <span dir="ltr"><<a href="mailto:pramsey@cleverelephant.ca" target="_blank">pramsey@cleverelephant.ca</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Move up to the latest patch revision of postgis, the 2.1.0 release<br>
definitely had a few subtle issues with geography distance of the sort<br>
you are describing.<br>
<br>
P.<br>
<div><div class="h5"><br>
On Thu, Jan 8, 2015 at 2:09 PM, Hartranft, Robert M.<br>
(GSFC-423.0)[RAYTHEON COMPANY] <<a href="mailto:robert.m.hartranft@nasa.gov">robert.m.hartranft@nasa.gov</a>> wrote:<br>
> Hi All,<br>
><br>
> I am hoping someone here has experience with the<br>
> issue I am having…<br>
><br>
><br>
> I am using:<br>
> -rw-r--r-- 1 postgres postgres 10241357 Oct 8 22:04 gdal-1.10.1.tar.gz<br>
> -rw-r--r-- 1 postgres postgres 1813726 Oct 8 22:04 geos-3.4.2.tar.bz2<br>
> -rw-r--r-- 1 postgres postgres 340953 Oct 8 22:04 json-c-0.9.tar.gz<br>
> -rw-r--r-- 1 postgres postgres 5161069 Oct 8 22:04 libxml2-2.9.0.tar.gz<br>
> -rw-r--r-- 1 postgres postgres 16930885 Oct 8 22:04 perl-5.16.3.tar.gz<br>
> -rw-r--r-- 1 postgres postgres 195442 Oct 8 22:04 pgtap-0.94.0.tar.gz<br>
> -rw-r--r-- 1 postgres postgres 6518378 Oct 8 22:04 postgis-2.1.0.tar.gz<br>
> -rw-r--r-- 1 postgres postgres 21865589 Oct 8 22:04<br>
> postgresql-9.3.4.tar.gz<br>
> -rw-r--r-- 1 postgres postgres 785279 Oct 8 22:04 proj-4.8.0.tar.gz<br>
><br>
> I have a table that contains a geography column.<br>
> The column contains “multipolygon” values.<br>
><br>
> When I do a query that uses the ST_Intersects() function it<br>
> seems to return the wrong results if the underlying plan uses<br>
> the _st_distance() function and it returns the correct results if<br>
> the plan uses the st_intersects() function.<br>
><br>
> I am having the issue with a multipolygon column but it might<br>
> happen of other types as well.<br>
><br>
><br>
><br>
> Here is an small example of the issue:<br>
><br>
> =# create table multipolygon<br>
> (<br>
> mpid bigint not null,<br>
> multip Geography(multipolygon, 4326) not null<br>
> );<br>
> CREATE TABLE<br>
><br>
><br>
> =# alter table multipolygon add constraint PK_POLYGON primary key (mpid);<br>
> ALTER TABLE<br>
><br>
><br>
> Now I will add two small “rectangles” into the table. I don’t really need<br>
> to<br>
> use multiple polygons to reproduce the issue so I will use a single polygon<br>
> with<br>
> points to approximate lines of constant latitude (just trying to avoid the<br>
> great<br>
> circle arc).<br>
><br>
> =# insert into multipolygon<br>
> values (1, ST_GeographyFromText('SRID=4326;multipolygon(((20 10, 21 10, 22<br>
> 10, 23 10, 24 10, 25 10, 26 10, 27 10, 28 10, 29 10, 30 10, 30 20, 29 20, 28<br>
> 20, 27 20, 26 20, 25 20, 24 20, 23 20, 22 20, 21 20, 20 20, 20 10)))'));<br>
> INSERT 0 1<br>
><br>
> =# insert into multipolygon<br>
> values (2, ST_GeographyFromText('SRID=4326;multipolygon(((25 15, 26 15, 27<br>
> 15, 28 15, 29 15, 30 15, 31 15, 32 15, 33 15, 34 15, 35 15, 35 25, 34 25, 33<br>
> 25, 32 25, 31 25, 30 25, 29 25, 28 25, 27 25, 26 25, 25 25, 25 15)))'));<br>
> INSERT 0 1<br>
><br>
><br>
><br>
> Now I run a query that should return both rows.<br>
><br>
><br>
> =# select m1.mpid from multipolygon m1 where st_intersects(m1.multip<br>
> ,'POLYGON((27 17, 28 17, 28 19, 27 19, 27 17))'::geography);<br>
> mpid<br>
> ------<br>
> 1<br>
> (1 row)<br>
><br>
> =# select m1.mpid from multipolygon m1 where st_intersects((select multip<br>
> from multipolygon m2 where m2.mpid = m1.mpid ),'POLYGON((27 17, 28 17, 28<br>
> 19, 27 19, 27 17))'::geography );<br>
> mpid<br>
> ------<br>
> 1<br>
> 2<br>
> (2 rows)<br>
><br>
><br>
> You can see the two queries seem to be logically the same but they return<br>
> different results<br>
> Next I will look at the plans for the two queries<br>
><br>
><br>
> =# explain select m1.mpid from multipolygon m1 where st_intersects(m1.multip<br>
> ,'POLYGON((27 17, 28 17, 28 19, 27 19, 27 17))'::geography);<br>
><br>
> QUERY PLAN<br>
> -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------<br>
> Seq Scan on multipolygon m1 (cost=0.00..317.40 rows=77 width=8)<br>
> Filter: ((multip &&<br>
> '0103000020E610000001000000050000000000000000003B4000000000000031400000000000003C4000000000000031400000000000003C4000000000000033400000000000003B4000000000000033400000000000003B400000000000003140'::geography)<br>
> AND (_st_distance(multip,<br>
> '0103000020E610000001000000050000000000000000003B4000000000000031400000000000003C4000000000000031400000000000003C4000000000000033400000000000003B4000000000000033400000000000003B400000000000003140'::geography,<br>
> 0::double precision, false) < 1e-05::double precision))<br>
> (2 rows)<br>
><br>
><br>
> =# explain select m1.mpid from multipolygon m1 where st_intersects((select<br>
> multip from multipolygon m2 where m2.mpid = m1.mpid ),'POLYGON((27 17, 28<br>
> 17, 28 19, 27 19, 27 17))'::geography );<br>
><br>
> QUERY PLAN<br>
> -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------<br>
> Seq Scan on multipolygon m1 (cost=0.00..9788.80 rows=387 width=8)<br>
> Filter: st_intersects((SubPlan 1),<br>
> '0103000020E610000001000000050000000000000000003B4000000000000031400000000000003C4000000000000031400000000000003C4000000000000033400000000000003B4000000000000033400000000000003B400000000000003140'::geography)<br>
> SubPlan 1<br>
> -> Index Scan using pk_polygon on multipolygon m2 (cost=0.15..8.17<br>
> rows=1 width=32)<br>
> Index Cond: (mpid = m1.mpid)<br>
> (5 rows)<br>
><br>
> You can see that the query using _st_distance is the one that returns a bad<br>
> result.<br>
> However I can also show _st_distance() works if it is limited to finding<br>
> just one<br>
> result. Here is the same query that didn’t work, however this time I add a<br>
> constraint<br>
> to limit the result set to one row.<br>
><br>
><br>
> =# select m1.mpid from multipolygon m1 where st_intersects(m1.multip<br>
> ,'POLYGON((27 17, 28 17, 28 19, 27 19, 27 17))'::geography) and m1.mpid = 1;<br>
> mpid<br>
> ------<br>
> 1<br>
> (1 row)<br>
><br>
> =# select m1.mpid from multipolygon m1 where st_intersects(m1.multip<br>
> ,'POLYGON((27 17, 28 17, 28 19, 27 19, 27 17))'::geography) and m1.mpid = 2;<br>
> mpid<br>
> ------<br>
> 2<br>
> (1 row)<br>
><br>
> =# explain select m1.mpid from multipolygon m1 where st_intersects(m1.multip<br>
> ,'POLYGON((27 17, 28 17, 28 19, 27 19, 27 17))'::geography) and m1.mpid = 2;<br>
><br>
> QUERY PLAN<br>
> -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------<br>
> Index Scan using pk_polygon on multipolygon m1 (cost=0.15..8.43 rows=1<br>
> width=8)<br>
> Index Cond: (mpid = 2)<br>
> Filter: ((multip &&<br>
> '0103000020E610000001000000050000000000000000003B4000000000000031400000000000003C4000000000000031400000000000003C4000000000000033400000000000003B4000000000000033400000000000003B400000000000003140'::geography)<br>
> AND (_st_distance(multip,<br>
> '0103000020E610000001000000050000000000000000003B4000000000000031400000000000003C4000000000000031400000000000003C4000000000000033400000000000003B4000000000000033400000000000003B400000000000003140'::geography,<br>
> 0::double precision, false) < 1e-05::double precision))<br>
> (3 rows)<br>
><br>
><br>
> This might be the wrong conclusion but it appears that<br>
> _st_distance() doesn’t work if the result set contains more<br>
> than one row.<br>
><br>
><br>
><br>
> Thanks in advance,<br>
> Bob<br>
> -----------------------------------------------------------------<br>
> Robert Hartranft<br>
> Mission Operations and Services<br>
> Intelligence & Information System<br>
> Raytheon Company<br>
> 301-851-8197<br>
> <a href="mailto:robert_m_hartranft@raytheon.com">robert_m_hartranft@raytheon.com</a><br>
> <a href="mailto:robert.m.hartranft@nasa.gov">robert.m.hartranft@nasa.gov</a><br>
> 5700 Rivertech Court<br>
> Riverdale, MD 20737<br>
> -----------------------------------------------------------------<br>
><br>
><br>
</div></div>> _______________________________________________<br>
> postgis-users mailing list<br>
> <a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
> <a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a></blockquote></div><br></div>