[prev in list] [next in list] [prev in thread] [next in thread] 

List:       geos-devel
Subject:    Re: [geos-devel] polar geometry spanning 0-360 lon
From:       David Shean <dshean () gmail ! com>
Date:       2014-03-11 22:37:52
Message-ID: 83657635-D74D-4812-B234-0D65968307DD () gmail ! com
[Download RAW message or body]

On Mar 11, 2014, at 2:53 PM, Mike Toews <mwtoews@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@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@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/geos-devel

_______________________________________________
geos-devel mailing list
geos-devel@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/geos-devel


[prev in list] [next in list] [prev in thread] [next in thread]