9

I am using the ST_Distance function to calculate the distance between two geometries (a train station and a building) Since I know that all buildings and all train stations are in Chicago, which has an excellent/complete street grid, I would like to use Manhattan (or taxicab) distance.

The generic formula for this is the difference in X plus the difference in Y, so Abs(X1-X2)+Abs(Y1-Y2).

What PostgreSQL query would make this work?

Improve this question
4
  • 1
    A quick thought here: Your grid 'x' and 'y' aren't necessarily aligned with the co-ordinate system's 'x' and 'y'. So you may need to rotate the vector before extracting the components and calculating.
    – Craig Ringer
    Commented Aug 13, 2014 at 16:07
  • @CraigRinger I transform the coordinates to the locally prevailing projection, EPSG 3435, Illinois StatePlane East Feet. This is used by the City of Chicago for all of its GIS work. I've answered my own question with some validation using Google Maps walking distance calculation.
    – stevevance
    Commented Aug 13, 2014 at 16:13
  • 1
    Have you also considered extending your PostGIS data with the pgRouting module and using it's built-in functions? Apparently the A* method uses something similar.
    – RyanKDalton
    Commented Aug 13, 2014 at 16:20
  • @RyanDalton I have considered using pgRouting for another project of mine but the hassle of setting one up for this project isn't worth the more accurate results or the resource cost in calculating the route.
    – stevevance
    Commented Aug 13, 2014 at 16:43

2 Answers 2

Reset to default
6

I am answering my own question with a proposed query. 

select *, ABS(x_permit-x_station)+ABS(y_permit-y_station) as manhattan FROM (SELECT
longname AS NAME,
lines AS metadata,
T .slug,
ST_Distance (
    T .geom,
    ST_Transform (P .geometry, 3435)
) AS distance, 
ST_X(ST_Transform(p.geometry, 3435)) as x_permit,
ST_Y(ST_Transform(p.geometry, 3435)) as y_permit,
ST_X(t.geom) as x_station,
ST_Y(t.geom) as y_station
FROM
permits P,
stations_cta T
WHERE
P .permit_ = '100533644'
ORDER BY
distance
LIMIT 2) as foo

This results in the following with some columns snipped out:

Kedzie-Ravenswood   Brown Line  3738.52830193659    3796.29623843171
Addison-O'Hare  Blue Line   4105.37381385087    5790.20002649655

The first numbered column is the distance (in feet, because I'm using EPSG 3435) calculated by the ST_Distance PostGIS function, and the second numbered column is the result of the Manhattan distance formula.

I spot-checked the second result, getting the walking distance from Google Maps between the Addison Blue Line CTA station and the building at 3226 W Belle Plaine Ave (noted as '100533644' in the query). Google Maps outputted a 1.1 miles walking route while the Postgres result of 5,790 feet = 1.09 miles. The difference is acceptable for my purposes.

Improve this answer
2
  • 1
    this is great - you might have solved a problem that is otherwise dealt with in PGRouting 'driving distance' function... I will test this out for some problems we have at DPS...!
    – Inactivated Account
    Commented Aug 13, 2014 at 17:10
  • @mapBaker Thanks! pgRouting is too complex for my simple mathematics need.
    – stevevance
    Commented Aug 13, 2014 at 18:50
3

I think I also found a slightly more elegant solution that uses trigonometry and the built in ST_Azimuth function and encapsulated it into a nice function:

CREATE OR REPLACE FUNCTION JZ_TaxiCab(p1 geometry, p2 geometry)
RETURNS REAL AS $$
DECLARE 
    az REAL;
    h REAL;
BEGIN
    az := ST_Azimuth(p1, p2);
    /* Cast to geography to get result in meters */
    h := ST_Distance(p1::geography, p2::geography);
    /*   Note: we have to take abs() because the values returned by
         can be positive or negative. We really don't necessarily care
         about the reference point since it's going to be a right triangle.
    */
    RETURN h * abs(sin(az)) + h * abs(cos(az));
END;
$$  LANGUAGE plpgsql
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.