A PostGIS user asked me an interesting question this week, looking for the intersection of these two polygons:
POLYGON((170 50,170 72,-130 72,-130 50,170 50))
POLYGON((-170 68,-170 90,-141 90,-141 68,-170 68))
Now, from the coordinate sizes, you can see these shapes are described in longitude/latitude. So naturally you should consider using the PostGIS “geography” type, particularly since the shapes are quite big.
However, when you plug the shape into SQL and try out the calculation:
SELECT ST_AsText(
ST_Intersection(
ST_GeogFromText('POLYGON((170 50,170 72,-130 72,-130 50,170 50))'),
ST_GeogFromText('POLYGON((-170 68,-170 90,-141 90,-141 68,-170 68))')
));
You get a bad result:
ERROR: transform: couldn't project point (-170 90 0): tolerance condition error (-20)
Hmm. Transform? This sounds like reprojection. We’re working in geography though, what’s going on? Check the definition of the ST_Intersects(geography, geography) function.
CREATE OR REPLACE FUNCTION ST_Intersection(geography, geography)
RETURNS geography
AS 'SELECT
geography(
ST_Transform(
ST_Intersection(
ST_Transform(
geometry($1),
_ST_BestSRID($1, $2)),
ST_Transform(
geometry($2),
_ST_BestSRID($1, $2))
), 4326))'
LANGUAGE 'SQL' IMMUTABLE STRICT;
Crazy stuff. Geographies are being cast to geometry, transformed into a planar projection (magically selected by _ST_BestSRID) intersected, then transformed back into geographics and cast back to geography.
The problem is that _ST_BestSRID is failing to pick an appropriate projection to do the calculation in. Now, there are some shapes that are just too big to pick any projection for, but in this case the shapes are a good deal smaller than a hemisphere. However, they are inconveniently located:
One of them goes all the way up to the pole, but the other extends moderately far south. The “best SRID” routine evaluates first whether the shapes can fit into a single UTM zone (they can’t), if they are sufficiently far north to use a polar stereographic (only one is), and finally falls back to a world mercator projection. But, one of the shapes goes up to the pole (which mercator places at infinity), so in the end even the mercator projection fails.
What to do? In order to accurately derive a new shape in planar space that reflects the result on the sphere, we need a projection that represents great circles (used to connect points on the sphere) as straight lines (used to connect points on the plane). We need a Gnomic projection.
Fortunately the projection library underneath PostGIS — Proj4 — supports the gnomic projection. Now, we need a gnomic projection that is parameterized to optimally fit our geography shapes. The approximate center of the two shapes falls at (-160, 70) so we’ll use that as the center of the projection.
Now we insert a record into the PostGIS “spatial_ref_sys” table, giving our new projection an identifier number of 55000:
INSERT INTO spatial_ref_sys (srid, proj4text)
VALUES (55000, '+proj=gnom +lat_0=70 +lon_0=-160');
Now we can run the intersection again, this time in our specially chosen gnomic:
WITH shapes AS
(
SELECT
ST_GeogFromText('POLYGON((170 50,170 72,-130 72,-130 50,170 50))') AS p1,
ST_GeogFromText('POLYGON((-170 68,-170 90,-141 90,-141 68,-170 68))') AS p2
)
SELECT ST_AsText(
geography(ST_Transform(
ST_Intersection(
ST_Transform(geometry(p1), 55000),
ST_Transform(geometry(p2), 55000)
),4326)))
FROM shapes;
And this time we get a result:
POLYGON((-169.908 74.035,-141.155 73.424,-141 68,-170 68,-169.908 74.035))
As you can see eyeballing the globe above the resultant shape should have four corners (check), the southern two corners should be the southern two corners of the second input shape (check), and the northern two corners should be derived from the intersection of the northern edge of the first shape and the eastern/western edges of the second (check).
Astute readers of the geometry text will note that the northern edge of the first shape runs from (170 72) to (-130 72): why doesn’t the north edge of the derived shape fall on the 72d parallel? Because the shortest path between those points doesn’t run along the parallel, it’s a great circle arc that stretches to the north along its path.
This example has made it clear to me that the “best SRID” code needs a little upgrade, so that the logic first checks if UTM is acceptable, then switches into finding a best gnomic, and only breaks down and tries mercator if the shape is drastically too large to fit in a hemisphere.
“Best SRID” finding will always remain a hack, and direct math on the sphere is preferable, but for now we only have a few operations implemented on the sphere (area, length, distance) so we need some hacks to get our work done.
