With our latest release of PostGIS (release 2.3.3), projection datum shift files have been included across all the platforms. You will find them located in: bigsql/pg96/share/postgresql/contrib/postgis-2.3/proj Why do you care? Well, this means you can easily reproject your PostGIS data using the ST_Transform command. > ST_Transform — Return a new geometry with its coordinates transformed to a different spatial reference.

Exercise Prerequisites

Previous knowledge of projections and spatial references are assumed for this tutorial. To learn more about projections, I recommend these resources:

Previous blog posts to help you get setup with PostGIS:


What we will do in this exercise

Suppose we would like to find all the subway entrances that are withing 220 meters (about 2 football fields as the crow flies) to the Empire State Building (located at 40.748855, -73.985773).

A good projection for measuring distances in NYC is based on the Transverse Mercator projection: UTM Zone 19N (srid 32619).

Unfortunately, the data available by the MTA, Subway Entrances, has a spatial reference of WSGS:84 (srid 4326). This is a good choice if you want a map which shows the whole world, but not for measuring distances within New York City.

So.. we need to:

  • import the shapefile into our PostGIS database
  • create a copy of the table and reproject it
  • select all subway entrances located within 220 meters of the Empire State building converted from lat, long (40.748855, -73.985773) to it’s location in UTM Zone 19N.

Import Shapefile to PostGIS

I’ve made it easier for you by providing the data for you to download here: subway_entrances.zip

Unzip the shapefile, navigate to the location of your subway_entrances.shp file, and run the following command:

shp2pgsql -I -s 4326 subway_entrances.shp subway_entrances | psql -U postgres -d postgis_sample

Connect to your database to make sure it was successfully imported. If using Windows, you may need to include -U postgres with the psql command:

psql 
postgres=# \connect postgis_sample;

See that subway_entrances table is in the database:

postgis_sample=# \dt

Now look at the attributes in subway_entrances table:

\d+ subway_entrances

Column Type Modifiers Storage Stats target Description
gid integer not null default nextval(‘subwayentrancesgid_seq’::regclass) plain
objectid numeric main
url character varying(254) extended
name character varying(254) extended
line character varying(254) extended
geom geometry(Point,4326) main

Reproject subway_entrances from srid 4326 to 32619

Create a table subway_entrances_utm that is a duplicate of subway_entrances:

CREATE TABLE subway_entrances_utm AS SELECT * FROM subway_entrances;

Reproject the new table to UTM Zone 19N (srid 32619):

ALTER TABLE subway_entrances_utm
    ALTER COLUMN geom TYPE geometry(POINT,32619) 
    USING ST_Transform(geom,32619);

Select stations 220 meters from Empire State Building

Run the following query composed of these PostGIS commands:

  • ST_DWithin
  • ST_Transform
  • ST_SetSRID
  • ST_MakePoint

    SELECT name, line
    FROM subway_entrances_utm
    WHERE ST_DWithin(geom,ST_Transform(ST_SetSRID(ST_MakePoint(-73.985773,40.748855),4326),32619),220);
    

Your result should show 7 rows, including one having no value for the name field:

name               |     line
---------------------------------+---------------
 6th Ave & 34th St at SE corner  | B-D-F-M-N-Q-R
 6th Ave & 34th St at NE corner  | B-D-F-M-N-Q-R
 Broadway & 32nd St at NE corner | B-D-F-M-N-Q-R
 Broadway & 32nd St at NW corner | B-D-F-M-N-Q-R
 6th Ave & 35th St at SE corner  | B-D-F-M-N-Q-R
 6th Ave & 35th St at NE corner  | B-D-F-M-N-Q-R
                                 | B-D-F-M-N-Q-R

If you want, you can create a new table, subway_entrances_empire, that only contains these subway stations:

CREATE TABLE subway_entrances_empire AS SELECT *
        FROM subway_entrances_utm
        WHERE ST_DWithin(geom,ST_Transform(ST_SetSRID(ST_MakePoint(-73.985773,40.748855),4326),32619),220);