[GIS] Importing MultiPolygon into PostGIS

postgis

I have a MultiPolygon with 2 relatively simple polygons in it:
http://microflush.org/json/MultiPolygon.json

I've pulled out the 2 polygons from the above MultiPolygon below, just to make sure they both exist (same url above, but different json):

json/upper.json
json/lower.json

Any/all of the 3 above can be copied/pasted into a GeoJSON tester to view them: http://geojsonlint.com/

I've stored the original MultiPolygon in PostgreSQL as below:

$ ogr2ogr -f "PostgreSQL" PG:"dbname=weatherzones user=postgres" "MultiPolygon.json" -nln polys

Here's how it looks in Postgres after the import:

polytest=# \d polys;

                                       Table "public.polys"
    Column    |          Type           |                        Modifiers                        
--------------+-------------------------+---------------------------------------------------------
 ogc_fid      | integer                 | not null default nextval('polys_ogc_fid_seq'::regclass)
 wkb_geometry | geometry(Geometry,4326) | 
 warnings     | character varying       | 
Indexes:
    "polys_pk" PRIMARY KEY, btree (ogc_fid)
    "polys_geom_idx" gist (wkb_geometry)

When I try to query postgres using single points which lie in the lower polygon, I do not get rows returned. However, when I query using points which lie in the upper polygon, I get the row returned.

ie:

This point lies in the lower polygon, and this query returns no rows:

select warnings from polys 
where ST_Intersects(ST_PointFromText('POINT( -116.024551 38.485773 )', 4326), wkb_geometry);

This next point lies in the upper polygon, and this query returns a row, which is the MultiPolygon:

select warnings from polys 
where ST_Intersects(ST_PointFromText('POINT( -114.879913 39.249129 )', 4326), wkb_geometry);

Is there something wrong with my query, the way I imported the data, or perhaps the original MultiPolygon which can explain why a query with a point from the lower polygon won't return the MultiPolygon row? I've done this with other MultiPolygons and have not had an issue like this.

Best Answer

Your MultiPolygon is invalid. You can test it with ogrinfo:

ogrinfo -ro -dialect sqlite -sql "select IsValid(geometry) from
OGRGeoJSON" multipolygon.json
INFO: Open of `multipolygon.json'
      using driver `GeoJSON' successful.
GEOS warning: Self-intersection at or near point -115.907104 39.162109000000001

Layer name: SELECT
Geometry: None
Feature Count: 1
Layer SRS WKT:
(unknown)
isvalid(geometry): Integer (0.0)
OGRFeature(SELECT):0
  isvalid(geometry) (Integer) = 0

You can try to correct it in PostGIS with ST_MakeValid http://postgis.net/docs/ST_MakeValid.html. However, it is not simple

Step 1 Correct the faulty geometry:

create table valid as select st_makevalid(geometry) as geometry from multipolygon;

Step 2 Test:

select * from valid where ST_Intersects(ST_PointFromText('POINT( -116.024551 38.485773 )',4326), geometry);

Result: error

ERROR:  Relate Operation called with a LWGEOMCOLLECTION type.  This is unsupported.
HINT:  Change argument 2: 'GEOMETRYCOLLECTION(POLYGON((-115.000793 38.677307,-115.000793 38.499878,-115....'

Pity that ST_Intersects does not handle geometrycollections. So you must get rid of the collection. First part of the collection is the polygon so you can take just that.

Step 3 Keep just the polygon part

update valid set geometry=ST_GeometryN(geometry,1);

Step 4 Test again 

select * from valid where ST_Intersects(ST_PointFromText('POINT( -116.024551 38.485773 )',4326), geometry);

Result: Success!

Related Question