2

When I perform an intersection of a polygon I draw in OpenLayers and a postgis database layer, it seems that I am getting incorrect results.

Intersection works correctly on some layers. For instance, if I intersect a triangle with a layer of polygons that represent crop fields, I get the following:

enter image description here

The query my app generates to produce the above result is:

SELECT ST_AsText(ST_Intersection(%(geometries_0)s::geometry, %(geometry)s::geometry))

where geometries_0 is my triangle:

POLYGON((-104.84928345939991 40.518951354186285,-104.82319093011056 40.51953858115158,-104.83700967095314 40.50707521626648,-104.84928345939991 40.518951354186285))

and geometry is my layer of crop fields, as well-known text:

MULTIPOLYGON(((-104.841309611298 40.5075331998226,-104.84173356681 40.5069932245841,-104.842041204329 40.50640946683,-104.842224948796 40.5057962996657,-104.842280275816 40.5051688207073,-104.842205823049 40.5045424803865,-104.842003423773 40.5039327015263,-104.841678061729 40.5033544995574,-104.841237748411 40.502822112724,-104.840693325791 40.5023486513933,-104.840058199365 40.5019457751149,-104.839348008051 40.5016234053897,-104.838580239118 40.5013894812384,-104.837773797582 40.5012497635973,-104.836948540713 40.501207693373,-104.836124789073 40.5012643066572,-104.83532282616 40.5014182091969,-104.834562398965 40.5016656107496,-104.833862231727 40.5020004184754,-104.833239564888 40.5024143870601,-104.832709730574 40.5028973218633,-104.83228577506 40.5034373300773,-104.831978137541 40.5040211136997,-104.831794393074 40.5046342970926,-104.831739066055 40.5052617810515,-104.831813518821 40.5058881146554,-104.832015918097 40.5064978757385,-104.832341280141 40.5070760506105,-104.83278159346 40.5076084036796,-104.833326016079 40.5080818278834,-104.833961142505 40.5084846673069,-104.834671333819 40.5088070040565,-104.835439102753 40.5090409023397,-104.836245544289 40.5091806037522,-104.837070801158 40.5092226689759,-104.837894552799 40.5091660624086,-104.83869651571 40.509012177646,-104.839456942906 40.5087648031902,-104.840157110143 40.5084300292289,-104.840779776982 40.5080160977709,-104.841309611298 40.5075331998226)))

However, if I perform the same query with a different layer ("soils"), I get an empty result:

enter image description here

The query is the same:

SELECT ST_AsText(ST_Intersection(%(geometries_0)s::geometry, %(geometry)s::geometry))

with a polygon geometries_0 that should overlap:

POLYGON((-104.84627938530097 40.54511058649626,-104.83460641167578 40.545175808723876,-104.84070039055733 40.537283458057615,-104.84627938530097 40.54511058649626))

and a geometry layer representing soils, similar to the crop fields in the above query:

MULTIPOLYGON(((-104.939716 40.258166,-104.939775 40.258174,-104.939963 40.258159,-104.940159 40.258065,-104.940039 40.257671,-104.939917 40.25749,-104.939928 40.257419,-104.94003 40.257404,-104.940265 40.257641,-104.940632 40.257902,-104.940826 40.258061,-104.941051 40.258188,-104.941123 40.258235,-104.941205 40.258283,-104.941246 40.258275,-104.941287 40.258212,-104.941186 40.258094,-104.941186 40.258007,-104.941167 40.257921,-104.941105 40.257858,-104.941044 40.257786,-104.941045 40.257716,-104.941127 40.257676,-104.94122 40.257653,-104.94141 40.257731,-104.941559 40.257671,-104.941255 40.257181,-104.940857 40.256794,-104.940644 40.256478,-104.940319 40.255997,-104.940003 40.255728,-104.939676 40.255561,-104.939419 40.255544,-104.938895 40.255529,-104.938287 40.255512,-104.938046 40.255528,-104.937733 40.255549,-104.937322 40.255533,-104.937012 40.255577,-104.936947 40.255593,-104.936623 40.255774,-104.936581 40.255924,-104.93658 40.256042,-104.936661 40.256223,-104.93671 40.256436,-104.936842 40.256618,-104.937262 40.256753,-104.937662 40.256818,-104.937897 40.257,-104.938181 40.25745,-104.938374 40.257742,-104.938465 40.257931,-104.938782 40.258051,-104.939121 40.258092,-104.939439 40.258133,-104.939716 40.258166)))

I use the postgis ST_AsText function to convert database layers to well-known text, and I checked to make sure that all layers have the EPSG:4326 projection (using the Find_SRID function).

Why would one layer (crop fields) intersect correctly, and another (soils) not? I've tried the same query using geographies instead of geometries, with the same result.

3
  • Don't cross post: stackoverflow.com/questions/18005963/…
    – Mike T
    Aug 1, 2013 at 23:34
  • 3
    This question has been asked and already answered on Stack Overflow - it looks like a GIS question so I would have preferred to see it posted here rather than there.
    – PolyGeo
    Aug 2, 2013 at 0:00
  • I've rewarded this (+1) in spite of voting to close it because it is an excellently presented question and I really would like to see the asker bring more such questions to GIS Stack Exchange.
    – PolyGeo
    Aug 2, 2013 at 0:25

1 Answer 1

1

As this question is a cross-post, I have an answer there: https://stackoverflow.com/a/18006366/327026

1
  • 1
    (+1) it's a shame we can't easily migrate your answer here and flag the question there as a duplicate since here seems the more logical home and it was a very well phrased question.
    – PolyGeo
    Aug 2, 2013 at 0:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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