1

I have two WGS84 NAD83 shapefiles, one polygonal, the other containing lines:

Geometry: Polygon
Feature Count: 1294
Extent: (-79.761059, 44.991143) - (-57.105486, 62.582865)
Layer SRS WKT:
GEOGCS["GCS_North_American_1983",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS_1980",6378137.0,298.257222101]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]

Geometry: Line String
Feature Count: 427227
Extent: (-79.882212, 44.929534) - (-55.834430, 62.427856)
Layer SRS WKT:
GEOGCS["GCS_North_American_1983",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS_1980",6378137.0,298.257222101]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]

which I export to PostGIS 2 using the -G flag, to obtain geographies:

$ shp2pgsql -G -I -s 4326 polygons | psql ...
$ shp2pgsql -G -I -s 4326 lines | psql ...

However when I perform a query to obtain their intersection:

select * from polygons p, lines l
where p.geog && l.geog 
and st_intersects(p.geog, l.geog)

it runs for a very, very long time, and doesn't yield the same results as the one performed over the equivalent tables, but with geometries instead (i.e. shp2pgsql without the -G flag), which runs very rapidly. According to the docs, st_intersects should play well with both types, so why this difference?

  • Can you add the extent of the two shapefiles? – AndreJ Mar 12 '14 at 5:15
  • They both span the entire Quebec province, I updated my question with their infos. – cjauvin Mar 12 '14 at 13:55
  • Woups, I just realized that I assumed my shapefiles to be WGS84, whereas they seem to be NAD83! I guess that explains my problem, as geographies work only with the latter, right? Shouldn't both systems be more or less equivalent though? – cjauvin Mar 12 '14 at 14:04
2

The results may very legitimately be different.

Interpretations on a plane are not the same as interpretations on a sphere, and that's what you're seeing when you compare geometry to geography calculations.

1

From the documentation:

For geography, this function has a distance tolerance of about 0.00001 meters and uses the sphere rather than spheroid calculation.

Unless the CRS you are using for your geometries also uses a sphere, the result will always be different.

1

The Postgis documentation notes that SRID 4326 is mandatory for geography, and calculations will take longer because they are done on a sphere instead of a plane as would be done for projected coordinates.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.