PostGIS intersection troubleshooting

Question

I just started using PostGIS & Postgresql and everything is running smoothly for the most part. When I try to find which MULTIPOLYGON s a POINT lies in I'm getting stuck. I have two separate points that I am certain lie inside one and only one shape that is of MULTIPOLYGON data type in my database. They are not the same points and they are in different formats.

Example 1, I'm not sure what format it is but the query returns true value like I expected (note, I found this value by loading the dataset into QGIS and hovering over a point inside).

In the second example, I geocoded an address that lies inside the shape I'm looking at. However, a false value is returned as a result of the query.

I used shp2pgsql to load the data into my database directly from a shape file. The shape file's SRID is 4269 (I've tried running the queries below while specifying the SRID during the GeomFromText call but the results are the same).

What is the difference between the two coordinate sets? What do I need to do so that I can perform an intersection test using POINTS that use lat/lon values?

1.) SELECT ST_Intersects((select the_geom from wardstable where gid=37), ST_GeomFromText('POINT(1172539 1924462)'));

2.) SELECT ST_Intersects((select the_geom from wardstable where gid=37), ST_GeomFromText('POINT(-87.6547884 41.96367)'));

Thanks!

Answer 1

Both the Multipolygons and the Points dataset should be in the same projection (SRID) when performing a spatial operation like ST_Intersects. In your second example the point's coordinates are in lat/lon (4326). You should transform them to 4269 using ST_Transform:

SELECT ST_Intersects((select the_geom from wardstable where gid=37), ST_Transform(ST_GeomFromText('POINT(-87.6547884 41.96367)',4326),4269));

Edit: I missed the SRID parameter in ST_GeomFromText.

Answer 2

I tried your suggestion a few different ways. First I checked out the SRID in the geometry_columns table. It was set to -1 because I believe I didn't specify the SRID in the shp2pgsql tool. But I tried your suggestion anyway a few different ways, all of them returning a false value...

SELECT ST_Intersects((select the_geom from wardstable where gid=37), 
ST_Transform(ST_GeomFromText('POINT(-87.6547884 41.96367)',4269),4269));  

SELECT ST_Intersects((select ST_Transform(the_geom,4269) from wardstable where gid=37), 
ST_Transform(ST_GeomFromText('POINT(-87.6547884 41.96367)',4269),4269));  

** I got an error unless I specified the SRID twice for the ST_GeomFromText used in conjuction with ST_Transform

I tried using SRID = 4326 just to see if I got any different results and nothing changed.

So I checked the documentation you linked to. It mentioned I had to have PROJ installed, so I did the version and everything checked out.

"POSTGIS="1.4.0" GEOS="3.2.2-CAPI-1.6.2" PROJ="Rel. 4.7.1, 23 September 2009" USE_STATS"

The last thing I tried was creating two new tables (output below) and running the above queries again (with correct parameters). I'm still getting false values. Any ideas on where I might be taking a wrong turn?

Note: I also posted the OGRINFO just incase that is helpful.

postgres@ubuntu:$  shp2pgsql -c -D -s 4269 -i -I Wards.shp public.wardstabledos | psql -d gisraw
Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]
SET
BEGIN
NOTICE:  CREATE TABLE will create implicit sequence "wardstabledos_gid_seq" for serial column "wardstabledos.gid"
NOTICE:  CREATE TABLE / PRIMARY KEY will create implicit index "wardstabledos_pkey" for table "wardstabledos"
CREATE TABLE
                         addgeometrycolumn                         
-------------------------------------------------------------------
 public.wardstabledos.the_geom SRID:4269 TYPE:MULTIPOLYGON DIMS:2 
(1 row)

CREATE INDEX
COMMIT
postgres@ubuntu:$  shp2pgsql -c -D -s 4326 -i -I Wards.shp public.wardstabletres | psql -d gisraw
Shapefile type: Polygon
Postgis type: MULTIPOLYGON[2]
SET
BEGIN
NOTICE:  CREATE TABLE will create implicit sequence "wardstabletres_gid_seq" for serial column "wardstabletres.gid"
NOTICE:  CREATE TABLE / PRIMARY KEY will create implicit index "wardstabletres_pkey" for table "wardstabletres"
CREATE TABLE
                         addgeometrycolumn                          
--------------------------------------------------------------------
 public.wardstabletres.the_geom SRID:4326 TYPE:MULTIPOLYGON DIMS:2 
(1 row)

CREATE INDEX
COMMIT

postgres@ubuntu:$ ogrinfo -so Wards.shp Wards
INFO: Open of `Wards.shp'
      using driver `ESRI Shapefile' successful.

Layer name: Wards
Geometry: Polygon
Feature Count: 53
Extent: (1091130.772400, 1813891.890100) - (1205199.881775, 1951669.020100)
Layer SRS WKT:
PROJCS["NAD_1983_StatePlane_Illinois_East_FIPS_1201_Feet",
    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]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",984250.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",-88.33333333333333],
    PARAMETER["Scale_Factor",0.999975],
    PARAMETER["Latitude_Of_Origin",36.66666666666666],
    UNIT["Foot_US",0.3048006096012192]]
DATA_ADMIN: Real (19.8)
PERIMETER: Real (19.8)......