polygons crossing 0/360 longitude

[Shane Byrne]

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@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Carlos Ferrão]

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@...> 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@...
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Shane Byrne]

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@...> 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@...
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>
> _______________________________________________
> postgis-users mailing list
> postgis-users@...
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

--

______________________________________________________________
Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
______________________________________________________________
Email: shane@...              Phone: [phone]
Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
______________________________________________________________
Mail : USGS - Astrogeology Division,
        2255 North Gemini Drive, Flagstaff, AZ 86001, US.
______________________________________________________________
_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Carlos Ferrão]

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@...> 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@...> 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@...
> >> http://postgis.refractions.net/mailman/listinfo/postgis-users
> >>
> > _______________________________________________
> > postgis-users mailing list
> > postgis-users@...
> > http://postgis.refractions.net/mailman/listinfo/postgis-users
> >
>
> --
>
> ______________________________________________________________
> Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
> ______________________________________________________________
> Email: shane@...              Phone: [phone]
> Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
> ______________________________________________________________
> Mail : USGS - Astrogeology Division,
>         2255 North Gemini Drive, Flagstaff, AZ 86001, US.
> ______________________________________________________________
> _______________________________________________
> postgis-users mailing list
> postgis-users@...
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[pcreso]

There is no easy way in PostGIS to do this that I'm aware of.

I've looked at scripts to do this

The approach I've taken that works fine in my situation, where the polygons are
background layers (with longitudes from -180 to 180) for mapserver & QGIS, is
to select all the polygons which intersect a -180/0/-90/90 polygon, translate
them 360 degrees, then add them to the original table.

This gives me a table which has a longitudinal extent of over -180 to 360
(depending on the extent of the intersecting polygons), with the W hemisphere
stored twice. Some data redundancy but a seamless coverage.

On a Linux system, I have used symbolic links to a raster tile set, to achieve
a similar -180-360 coverage. I have world files which apply to the links, so in
this case I have not needed to replicate the imagery for the W hemisphere, as
the links do this for me.

I don't know if this approach is any use to others, but it works for me
(sitting in New Zealand where +-180 is a major inconvenience).


Cheers,

  Brent Wood

_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Shane Byrne]

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@...> 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@...> 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@...
>> >> http://postgis.refractions.net/mailman/listinfo/postgis-users
>> >>
>> > _______________________________________________
>> > postgis-users mailing list
>> > postgis-users@...
>> > http://postgis.refractions.net/mailman/listinfo/postgis-users
>> >
>>
>> --
>>
>> ______________________________________________________________
>> Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
>> ______________________________________________________________
>> Email: shane@...              Phone: [phone]
>> Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
>> ______________________________________________________________
>> Mail : USGS - Astrogeology Division,
>>         2255 North Gemini Drive, Flagstaff, AZ 86001, US.
>> ______________________________________________________________
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users@...
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>
> _______________________________________________
> postgis-users mailing list
> postgis-users@...
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

--

______________________________________________________________
Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
______________________________________________________________
Email: shane@...              Phone: [phone]
Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
______________________________________________________________
Mail : USGS - Astrogeology Division,
        2255 North Gemini Drive, Flagstaff, AZ 86001, US.
______________________________________________________________
_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Stephen Woodbridge]

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@...> 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@...> 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@...
>>> >> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>> >>
>>> > _______________________________________________
>>> > postgis-users mailing list
>>> > postgis-users@...
>>> > http://postgis.refractions.net/mailman/listinfo/postgis-users
>>> >
>>>
>>> --
>>>
>>> ______________________________________________________________
>>> Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
>>> ______________________________________________________________
>>> Email: shane@...              Phone: [phone]
>>> Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
>>> ______________________________________________________________
>>> Mail : USGS - Astrogeology Division,
>>>         2255 North Gemini Drive, Flagstaff, AZ 86001, US.
>>> ______________________________________________________________
>>> _______________________________________________
>>> postgis-users mailing list
>>> postgis-users@...
>>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>>
>> _______________________________________________
>> postgis-users mailing list
>> postgis-users@...
>> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>
>

_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Carlos Ferrão]

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@...> 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@...> 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@...> 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@...
> >> >> http://postgis.refractions.net/mailman/listinfo/postgis-users
> >> >>
> >> > _______________________________________________
> >> > postgis-users mailing list
> >> > postgis-users@...
> >> > http://postgis.refractions.net/mailman/listinfo/postgis-users
> >> >
> >>
> >> --
> >>
> >> ______________________________________________________________
> >> Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
> >> ______________________________________________________________
> >> Email: shane@...              Phone: [phone]
> >> Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
> >> ______________________________________________________________
> >> Mail : USGS - Astrogeology Division,
> >>         2255 North Gemini Drive, Flagstaff, AZ 86001, US.
> >> ______________________________________________________________
> >> _______________________________________________
> >> postgis-users mailing list
> >> postgis-users@...
> >> http://postgis.refractions.net/mailman/listinfo/postgis-users
> >>
> > _______________________________________________
> > postgis-users mailing list
> > postgis-users@...
> > http://postgis.refractions.net/mailman/listinfo/postgis-users
> >
>
> --
>
> ______________________________________________________________
> Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
> ______________________________________________________________
> Email: shane@...              Phone: [phone]
> Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
> ______________________________________________________________
> Mail : USGS - Astrogeology Division,
>         2255 North Gemini Drive, Flagstaff, AZ 86001, US.
> ______________________________________________________________
> _______________________________________________
> postgis-users mailing list
> postgis-users@...
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>

_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[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@...> 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@...> 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@...
http://postgis.refractions.net/mailman/listinfo/postgis-users

[Shane Byrne]

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@...> 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@...> 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@...> 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@...
>>> >> >> http://postgis.refractions.net/mailman/listinfo/postgis-users
>>> >> >>
>
>

--

______________________________________________________________
Shane Byrne  -  University of Arizona  - Lunar & Planetary Lab
______________________________________________________________
Email: shane@...              Phone: [phone]
Web  : http://www.gps.caltech.edu/~shane  Fax  : [phone]
______________________________________________________________
Mail : USGS - Astrogeology Division,
       2255 North Gemini Drive, Flagstaff, AZ 86001, US.
______________________________________________________________

_______________________________________________
postgis-users mailing list
postgis-users@...
http://postgis.refractions.net/mailman/listinfo/postgis-users