1

I am very new to PostGIS (started using it this morning).

My goal is to determine if a lat/long point is contained inside a polygon that is described in a shapefile (The shapefile has a single polygon).

My plan for attack is to import the shapefile into PostGIS (this is the tool I have to use) and run a spatial query. These are my steps:

1) Infer the SRID from the shapefile:

cat CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.prj PROJCS["Calgary_3TM_WGS_1984_W114",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["World_Geodetic_System_of_1984_GEM_10C",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-114],PARAMETER["scale_factor",0.9999],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]

and http://spatialreference.org/ref/epsg/3776/

My best guess is that the SRID I should use is 3776. I have no background in GIS so I am just making a best guess.

2) Convert the shapefile into .sql file:

shp2pgsql -g the_geom_3776 -s "3776" -W "latin1" CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.shp > CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.sql

3) import the shapefile in PostGIS.

psql -d testdb -f CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.sql

I want to determine if a lat/long point is inside the perimeter of the City of Calgary. This are the coordinates: 51.0644894,-114.081923. Google tells me they are within the city limits so it must be true.

However, the following query returns false:

SELECT ST_Contains((select the_geom_3776 from calgis_citybound_limit), ST_SetSRID(ST_MakePoint(51.0644894,-114.081923),   3776));
st_contains 
-------------
 f
(1 row)

I have a gut feeling that my projection is wrong but I can't figure out what the problem might be.

  • The projection of your point is 4326 – Ian Turton Mar 6 '16 at 13:07
  • Thanks @iant. I think what you are suggesting is: SELECT ST_Contains((select the_geom_3776 from calgis_citybound_limit), ST_SetSRID(ST_MakePoint(51.0644894,-114.081923), 4326)); st_contains ------------- f (1 row) which doesn't seem to work. – nayden Mar 6 '16 at 14:10
3

Your coordinates are in lat/lon (EPSG:4326) but you need to convert them into the projection that your polygons are stored in, so your query should be something like:

SELECT ST_Contains((select the_geom_3776 from calgis_citybound_limit), 
    ST_Transform(
       ST_SetSRID(ST_MakePoint(51.0644894,-114.081923),4326),
    3776));
  • Thanks @iant. This did it. I was getting an error: ERROR: transform: couldn't project point (51.0645 -114.082 0): latitude or longitude exceeded limits (-14) but was able to fix it by switching the long/lat coordinates around. – nayden Mar 6 '16 at 14:25
1

Here are a few fixes to the query:

  1. Always use lng/lat (X/Y) axis order for PostGIS.
  2. As @iant mentioned, use ST_Transform to do the spatial transformation to 3TM.
  3. Flatten the query by removing the subquery (i.e. select (select ...)). This subquery is unnecessary and can be slower, as it bypasses any spatial index (if it exists).

This should work:

SELECT ST_Contains(the_geom_3776,
    ST_Transform(ST_SetSRID(ST_MakePoint(-114.081923, 51.0644894), 3776), 4326))
FROM calgis_citybound_limit;

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.