[postgis-users] polygons crossing 0/360 longitude

[Shane Byrne]

Hi Matt
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 
works.

Shane


CREATE OR REPLACE FUNCTION lon_wrap(GEOMETRY) RETURNS GEOMETRY AS $$
DECLARE
  end_geom GEOMETRY;
  line_geom GEOMETRY;
  xp FLOAT;
  yp FLOAT;
  nn INTEGER;

BEGIN
  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;
   END IF;
   SELECT INTO line_geom setpoint(line_geom,nn-1,MakePoint(xp,yp));
  END LOOP;

  SELECT INTO end_geom multi(MakePolygon(line_geom));

  SELECT INTO end_geom geomunion(
            multi(intersection( multi(box2d('BOX(-180 -90, 180 
90)')),end_geom )),
  translate(multi(intersection( multi(box2d('BOX(180 -90, 720 
90)')),end_geom )),-360,0,0)
  );

  ELSE
   end_geom := $1;
  END IF;

  RETURN end_geom;

END;
$$ 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,
> 
> Matt
> 
> ------------------------------------------- 
> NERC Earth Observation Data Centre 
> Space Science & Technology Dept. 
> CCLRC Rutherford Appleton Laboratory 
> Chilton, Didcot OX11 0QX 
> 01235 445645 (tel) / 445848 (fax) 
> http://www.neodc.rl.ac.uk/ 
> 

-- 

______________________________________________________________
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.
______________________________________________________________