How do I calculate Manhattan distance with PostGIS?
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.
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