[postgis-users] polygons crossing 0/360 longitude
Shane Byrne
shane at quake.mit.edu
Fri Jun 23 08:22:17 PDT 2006
Thanks a lot to all who responded. I've come up with a solution to my
problem. I wrote a plpgsql script to split polygons that cross the prime
meridian into multipolygons. It will handle complex polygons that cross
this meridian several times. I did not add the loop to allow it to
handle a multipolygon input (wasn't necessary for my case) although that
shouldn't be too hard. The code is copied below in case it can save
someone some time in the future.
Cheers
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 < 180) 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( 0 -90, 360
90)')),end_geom )),
translate(multi(intersection( multi(box2d('BOX(360 -90, 720
90)')),end_geom )),-360,0,0)
);
ELSE
end_geom := $1;
END IF;
RETURN end_geom;
END;
$$ LANGUAGE 'plpgsql';
>
> Carlos Ferrão wrote:
>> Hello Shane,
>>
>> One thing I do is to keep the original polygon coordinates to show to
>> the user. The longitude values are, in principle, limited between -180
>> and 180 on the interface. The search will work good and always
>> retrieve only one record even if internally you have two polygons
>> which intersect if someday -180 connects to 180.My approach is to be
>> pragmatic and deliver reasonably good solutions on time. Of course, if
>> the requirements change, this is no longer an option.
>>
>> The problem with your approach is that you're applying it to my
>> trivial example where, conveniently, the latitude is the same. Now
>> imagine that instead of 170 20, -170 20 you have 170 20, -170(=190)
>> 35. Imagine also that a polygon can have lines that cross the dateline
>> many times. Well, for each segment of those, you need to calculate m
>> (the declive) of the equation of the line y=mx+b and determine what
>> the latitude value will be when the longitude is 180. And this is if
>> you're plotting the polygons with lines, if you're using archs it gets
>> much more complicated.I would classify this as a last hopeless
>> solution.
>>
>> If you really want to keep your polygons neat, it's easier to use
>> geometry operations but you need to calculate the two polygons as I
>> said in my previous email, intersect them with some rectangle on the
>> <-180 >180 limits and get the complementary geometry.Something like
>> that should work.
>>
>> I use PERL for two reasons. The first is that I need to parse a file
>> with metadata of satellite products and PERL is very powerful with
>> regular expressions and very fast working with strings. The other
>> reason is that I prefer to process my data before I load it and I
>> don't want to put the responsibility of calculating coordinates inside
>> a specific database engine.
>>
>> Good luck,
>> Carlos.
>>
>> On 6/21/06, Shane Byrne <shane at quake.mit.edu> wrote:
>>> ok
>>> I was thinking more of that polygon (170 20, -170 20, 175 5, 170 10,
>>> 170
>>> 20) becoming:
>>>
>>> Polygon 1:
>>> 170 20, 180 20, 180 10, 175 5, 170 10, 170 20
>>> ^^^^^^ ^^^^^^
>>> Polygon 2:
>>> -180 20, -170 20, -180 10
>>> ^^^^^^^ ^^^^^^^
>>>
>>> i.e. introducing two new vertices on the dividing meridian and having
>>> the area of the two polygons add up to the total area of the original.
>>>
>>> It sounds like there is no easy way to do this in PostGis though. I
>>> think your approach of using perl to figure it out would probably be
>>> the
>>> easiest way.
>>>
>>> My polygons are pretty simple and only cross the meridian once, I guess
>>> a complicated polygon could do that several times which will confuse
>>> things.
>>>
>>> Cheers,
>>> Shane
>>>
>>>
>>> Carlos Ferrão wrote:
>>> > To calculate the two polygons I am using a PERL script on data
>>> > ingestion to build the insert statement..
>>> > If I have a polygon which is (lon/lat)
>>> > 170 20, -170 20, 175 5, 170 10, 170 20
>>> > it becomes
>>> > polygon 1:
>>> > 170 20, (-170+360=190) 20, 175 5, 170 10, 170 20
>>> >
>>> > polygon 2:
>>> > (170-360=-190) 20, -170 20, (175-360=-185) 5, (170-360=-190) 10,
>>> > (170-360=-190) 20
>>> >
>>> > Just insert the two polygons into a single row using the MULTIPOLYGON
>>> > object.
>>> >
>>> > Carlos.
>>> >
>>> > On 6/21/06, Shane Byrne <shane at quake.mit.edu> wrote:
>>> >> Thanks Carlos, this is a good idea.
>>> >>
>>> >> Finding the problem polygons is easy enough, but how did you
>>> split them
>>> >> into two polygons along this meridian? Can I define a line and
>>> with that
>>> >> use some postgis function to split the polygon?
>>> >> I'd also like to copy over all the information in the other
>>> fields to
>>> >> each new polygon.
>>> >>
>>> >>
>>> >> Shane
>>> >>
>>> >> Carlos Ferrão wrote:
>>> >> > Hi,
>>> >> > I recently had a similar problem with products crossing the
>>> -180/180
>>> >> > longitude. The problem is that postgis doesn't connect the -180
>>> to the
>>> >> > 180 as, for instance, Oracle Spatial Module.
>>> >> > The way I found to solve it was to do an algorithm to find the
>>> >> > polygons that cross the date line. For each one, I calculate 2
>>> >> > polygons, one in the negative coordinates and another in the
>>> positive
>>> >> > coordinates (you need to add/subtract 360). I add a
>>> multipolygon with
>>> >> > the information of the two polygons in a single row in the
>>> database
>>> >> > and it works perfectly.
>>> >> >
>>> >> > Hope it helps,
>>> >> > Carlos Ferrao
>>> >> > EOP - ESRIN - European Space Agency
>>> >> > Critical Software - www.criticalsoftware.com
>>> >> >
>>> >> > On 6/21/06, Shane Byrne <shane at quake.mit.edu> wrote:
>>> >> >> Hi,
>>> >> >> I have a postgres/postgis database with many (~80k) polygons, the
>>> >> >> vertices of which are stored in longitude (0-360) and latitude
>>> space.
>>> >> >> These polygons are typically small, but a few hundred of them
>>> span the
>>> >> >> prime meridian. This creates a problem for me when I search for
>>> >> polygon
>>> >> >> intersections. i.e. a polygon with a longitude range -5 to +5 is
>>> >> >> actually stored as +355 to +5 and so postgis thinks it
>>> intersects a
>>> >> >> whole bunch of polygons that it really doesn't. I could
>>> change to the
>>> >> >> -180 to 180 longitude system but that only moves the problem to a
>>> >> >> different longitude.
>>> >> >>
>>> >> >> Is there a graceful way around this?
>>> >> >> I'm thinking of just making two tables, the other table would
>>> have
>>> >> >> longitudes ranging from -180 to +180 and doing two searches with
>>> >> >> longitude ranges chosen to avoid the problem area on each
>>> table. Does
>>> >> >> postgis have a way to handle cyclical coordinates?
>>> >> >>
>>> >> >> Thanks for any help,
>>> >> >> Shane
>>> >> >>
>>> >> >>
>>> >> >> _______________________________________________
>>> >> >> postgis-users mailing list
>>> >> >> postgis-users at postgis.refractions.net
>>> >> >> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>> >> >>
>
>
--
______________________________________________________________
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
mailing list