[postgis-users] polygons crossing 0/360 longitude
shane at quake.mit.edu
Fri Aug 4 09:31:14 PDT 2006
I'm glad this has been useful to you. I think your problem is not caused
by how you are searching, but instead by how you are splitting the polygons.
The code I posted was suitable for polygons that span the 0/360 meridian
but won't work with polygons that span the -180/180 meridian. Here is a
version with a few numbers changed to work with your case. I haven't
been able to test this after making these changes so let me know if it
CREATE OR REPLACE FUNCTION lon_wrap(GEOMETRY) RETURNS GEOMETRY AS $$
IF ((xmax($1)-xmin($1)) > 180) THEN
SELECT INTO line_geom geometryn(boundary($1),1);
FOR nn IN 1 .. npoints(line_geom) LOOP
SELECT INTO xp x(pointn(line_geom,nn));
SELECT INTO yp y(pointn(line_geom,nn));
IF (xp < 0) THEN
xp := xp + 360;
SELECT INTO line_geom setpoint(line_geom,nn-1,MakePoint(xp,yp));
SELECT INTO end_geom multi(MakePolygon(line_geom));
SELECT INTO end_geom geomunion(
multi(intersection( multi(box2d('BOX(-180 -90, 180
translate(multi(intersection( multi(box2d('BOX(180 -90, 720
end_geom := $1;
$$ LANGUAGE 'plpgsql';
Pritchard, MJ (Matt) wrote:
> Hi Shane (& Carlos),
> Thanks for your posts about this: this is really helpful as I have been
> grappling with this issue, too. However I'm still stuck...
> Do you have an example of a spatial search query that works, Shane,
> using geometries loaded using your lon_wrap() function?
> The problem I have had is false hits caused by postgis incorrectly
> assigning a bounding box to a geometry spanning the dateline, which then
> spans the *rest* of the longitude range. E.g. if your geometry's
> longitude range is 170 to -170, its resulting BB will be -170 to 170
> (implying it crosses the Greenwich meridian, not the dateline), hence
> the false positive if you search over Greenwich.
> E.g. If you take a geometry like
> 'POLYGON(((-170 20,170 20,170 5,-170 5,-170 20)),-1'
> (which spans the dateline)
> And split it as per your method into
> 'MULTIPOLYGON(((360 20,360 5,190 5,190 20,360 20)),((0 5,0 20,170 20,170
> 5,0 5))),-1'
> Then apply a search like
> select * from matt_test where geom1 && geometryfromtext('POLYGON(( -10
> 10, 10 10, 10 5, -10 5, -10 10))',-1);
> then I still get a false positive. AFAICT, the problem is present
> whether your geometries are all in SRID=4326 or all in SRID=-1
> Have you got some other method of searching that avoids this problem?
> Many thanks,
> NERC Earth Observation Data Centre
> Space Science & Technology Dept.
> CCLRC Rutherford Appleton Laboratory
> Chilton, Didcot OX11 0QX
> 01235 445645 (tel) / 445848 (fax)
Shane Byrne - University of Arizona - Lunar & Planetary Lab
Email: shane at lpl.arizona.edu Phone: (928)556-7235
Web : http://www.gps.caltech.edu/~shane Fax : (928)556-7014
Mail : USGS - Astrogeology Division,
2255 North Gemini Drive, Flagstaff, AZ 86001, US.
More information about the postgis-users