[geos-devel] polar geometry spanning 0-360 lon

[David Shean]

On Mar 11, 2014, at 2:53 PM, Mike Toews <mwtoews at gmail.com> wrote:

> Hi David,
> 
> GEOS (and JTS) use Cartesian space, so concepts of latitude and
> longitude are irrelevant.

Ah, good to know.

> The choice of projection and precision are
> important for spatial operations, but these need to be handled
> externally using libraries like PROJ.4, or other logic for handling
> latitude and longitude.
> 
> The two valid polygons look like skinny lines in a Cartesian
> representation of EPSG:4326, but look more like squares when
> transformed to EPSG:3031. It's not clear what you were expecting, so I
> cant tell what is considered a "fail". The intersection results are
> the same for JTS, and does not seem like a bug to me.
> 
> Usually the best approach is to do spatial operations with GEOS in a
> projected Cartesian space (on either pole, for your case). But only do
> this if it is within the range of the projection's limits. This is
> further limited by the size of the geometry, as sometimes they don't
> fit within the projection nicely. But these are general GIS concerns
> that are out of scope for GEOS to consider.

This all makes sense.  Early on, I implemented the conversion to wgs84 for intersection operations to handle near-global geometries that exceeded the limits of the polar stereographic projections.  It sounds like I should always perform the intersection in projected coords, and handle these cases separately.  

Thanks for the quick response.

> 
> -Mike
> 
> On 12 March 2014 08:29, David Shean <dshean at gmail.com> wrote:
>> Hi list,
>> I've run into this a few times now, and I'm hoping someone can point me in the right direction.  I primarily use GEOS through the ogr API to work with polar data in either EPSG:3031 or EPSG:3413.  When I work with smaller geometries in these projections, I always convert to WGS84 before performing any operations, and this works well.  However, even this seems to fail for large polygons that span 0-360 around the pole:
>> 
>> g1 = POLYGON ((-45.0 -54.673528789190144,45.00716107726717 -54.669377943743747,135.0 -54.66522770972778,-135.00716107726717 -54.669377943743747,-45.0 -54.673528789190144))
>> g2 = POLYGON ((-45.0 -54.673528789190144,45.00255773311347 -54.672046393296853,135.0 -54.670564075397799,-135.002557733113463 -54.672046393296853,-45.0 -54.673528789190144))
>> g1.Intersection(g2)
>> POLYGON ((3.650402839379899 -54.671285186562592,-45.0 -54.673528789190144,-83.299396823994513 -54.67176254220486,3.650402839379899 -54.671285186562592))
>> 
>> For reference, the EPSG:3031 geom is POLYGON ((-3333500 3333500,3333500 3333500,3333500 -3333500,-3333500 -3333500,-3333500 3333500)).
>> 
>> This looks like a bug to me, but maybe there is a way to specify that the polygon should include minlat -90?  Thanks.
>> -David
> _______________________________________________
> geos-devel mailing list
> geos-devel at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/geos-devel