Sunday, 1 October 2017

postgis - Geometry vs geography for line-to-polygon intersection


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?




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...