0

I have two shapefiles and i'm trying to do simple spatial analysis. When I want to know which polygon belongs to other polygon in other shapefile I always use select by location and in Spatial selection method i choose contain the source layer feature.

Now i have these two shapefiles in postgis and i'm trying to do the same using SQL Query but i got these results:

Select * from state where ST_Contains(geom, (select geom from region where REGION_N='CAPITAL_REGION'));

=0 results, in ArcMap i got 1 result.

I also tried

ST_ContainsProperly

=0 results

When i tried

ST_Intersects

= 3 results, the same in ArcMap. So it looks like intersects works but "ST_Contains" doesn't work for some reason.

Where could be with SQL problem?

  • Could you post as WKT 1) REGION_N='CAPITAL_REGION' and 2) the one result polygon you get from ArcMap? – user30184 Jul 25 '15 at 21:19
  • @user30184 I used "ST_AsEWKT" with these results: 1) ERROR: WKB structure does not match expected size! CONTEXT: SQL function "st_astext" statement 1 the same result i got using "ST_AsText" 2) The result polygon from ArcMap "MULTIPOLYGON (((103.82894961825984 1.240635969100187,103.82898942561337 1.240622224623671,103.82902854737875...... 1.240641704556595,103.82894961825984 1.240635969100187)),((103.82233409716753 1.247288081083784,103.82235533213685 1.2,1,STATE_1 " Can't write whole result here since it's too long. – L.s. Jul 26 '15 at 10:30
  • Anyway i can provide WKT REGION_N='CAPITAL_REGION' from ArcMap too. I'm just not sure why i didn't get this result from PostGIS. I got error message you can see above. So here is the result: POLYGON ((103.84924291873233 1.362752820720889,103.84935645049903 1.362682483934137,103.84973091592528 1.362406260504348,103.84992386166969...... 1.353352855081298,103.82898956487392 1.353407496501254,103.82883357311633 1.353461875646926,103.82867779071746,1,CAPITAL_REGION – L.s. Jul 26 '15 at 11:00
  • Sorry, without complete geometries I can't guess what happens. How did you import your shapefiles into PostGIS? – user30184 Jul 26 '15 at 12:43
  • @user30184 Unfortunately it's not possible to add complete geometries here. Maybe could i send you complete geometries via message somehow? I used shp2pgsql. The exact command i used look like this: shp2pgsql -I -s 4326 PATH_TO\state.shp state | psql -U postgres -d MyDatabase. – L.s. Jul 26 '15 at 13:30
2

I could reproduce thew issue with your shapefile data. I used ogr2ogr for converting the shapefiles into PostGIS tables and your query finds nothing.

Select * from state where ST_Contains(geom, (select geom from region where REGION_N='CAPITAL_REGION'));

Next I opened your shapefiles with OpenJUMP and made some tests. Both layers have some topology errors (ring self-intersections) which can make trouble for spatial functions but in this case I think that the issue comes from very little differences in the coordinates of corresponding vertices in those two shapefiles. Here is an example:

Layer: region
Feature ID: 50
[      0:44] POINT (103.86042442121341 1.343437997354843)
[      0:45] POINT (103.86043696024463 1.3434034315038)

Layer: State
Feature ID: 48
[    7:0:45] POINT (103.86042442121337 1.3434379973548598)
[    7:0:46] POINT (103.86043696024463 1.3434034315038228)

The coordinates differ in 12th coordinate place or so which is very very little but they still differ. When PostGIS is making computations I believe that it finds that some vertices of CAPITAL_REGION are outside the polygon of STATE_1 and therefore ST_Contains if "false".

What I do usually in cases like this is to reduce the accuracy of vertices first to some reasonable measure. Six or seven decimals could be a good start (see https://en.wikipedia.org/wiki/Decimal_degrees).

The main reason for the trouble is perhaps that vertices which have the same coordinates in ArcMap get written with slightly different values into "State" and "regions" shapefiles. Or perhaps ArcMap is handling internally coordinates with some tolerance but PostGIS is using the full accuracy in analysis.

  • Thank you so much for this detailed explanation, I appreciate it. Yeah that's probably the main reason why ArcMap found 1 result in this case. I've never met with something like this, reduce the accuracy of vertices could be definitely solution. – L.s. Jul 26 '15 at 19:10

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.