[postgis-users] polygons crossing 0/360 longitude

[Shane Byrne]

ok, I think I'm half-way there...

Consider the polygon:
POLYGON((350.0 80.0,  10.0 80.0,  10.0 70.0, 355.0 70.0, 355.0 60.0,  
10.0 60.0,  10.0 50.0, 350.0 50.0, 350.0 80.0))

It's more complicated than my usual cases in that is crosses the prime 
meridian twice .

Fist unwrap the polygon so that the longitudes are continuous... add 360 
to x-vertices less than 180 [this is the step I need help with].
POLYGON((350.0 80.0, 370.0 80.0, 370.0 70.0, 355.0 70.0, 355.0 60.0, 
370.0 60.0, 370.0 50.0, 350.0 50.0, 350.0 80.0))

The following command would take that polygon and, as Stephen suggested, 
intersect it with boxes on either side of the prime meridian. then we 
translate the 360-720 longitude group back to the 0-360 range and then 
collect the results into a single multipolygon with geomunion.

select Astext(geomunion(
          multi(intersection(
GeomFromText('MULTIPOLYGON(((0.0 90.0,0.0 -90.0,360.0 -90.0, 360.0 
90.0,0.0 90.0)))'),
GeomFromText('MULTIPOLYGON(((350.0 80.0, 370.0 80.0, 370.0 70.0, 355.0 
70.0, 355.0 60.0, 370.0 60.0, 370.0 50.0, 350.0 50.0, 350.0 80.0)))')
)),
translate(multi(intersection(
GeomFromText('MULTIPOLYGON(((360.0 90.0,360.0 -90.0,720.0 -90.0, 720.0 
90.0,360.0 90.0)))'),
GeomFromText('MULTIPOLYGON(((350.0 80.0, 370.0 80.0, 370.0 70.0, 355.0 
70.0, 355.0 60.0, 370.0 60.0, 370.0 50.0, 350.0 50.0, 350.0 80.0)))')
)),-360,0,0)
));


This will produce a three element multipolygon with the proper vertices.

MULTIPOLYGON(((360 60,360 50,350 50,350 80,360 80,360 70,355 70,355 
60,360 60)),((0 70,0 80,10 80,10 70,0 70)),((0 50,0 60,10 60,10 50,0 50)))

I need help with the first step i.e. adding 360 to x-vertices less than 
180 to unwrap the polygon.
I'm a postgres novice [if that isn't already painfully obvious :) ].

Thanks,
Shane



Stephen Woodbridge wrote:
> Shane,
>
> If you take the polygon bbox and split that into two bboxes on 
> opposite sides of the -180/180 line, then you should be able to 
> intersect the original polygon with each of the right/left bbox 
> polygons to split it into the respective right/left components.
>
> -Steve
>
> Shane Byrne 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
>>>> >>