I have Point geometries created from lat,lon I loaded all the 10,000 stops in a certain city using gtfs

i would like to query my db to find the closest stop reference to my point.

the staypoint table is

    CREATE TABLE stay_points
        ad_id uuid NOT NULL,
        latitude numeric NOT NULL,
        longitude numeric NOT NULL,
        start_utc_timestamp timestamp without time zone NOT NULL,
        end_utc_timestamp timestamp without time zone NOT NULL,
        cluster_id numeric NOT NULL,
        radius numeric NOT NULL,
        geom geometry(Point,4326)
    WITH (
        OIDS = FALSE
    TABLESPACE pg_default;

and the stop table is :

    CREATE TABLE public.stops
    stop_name text COLLATE pg_catalog."default",
    stop_lon real NOT NULL,
    stop_lat real NOT NULL,
    stop_desc text COLLATE pg_catalog."default",
    stop_code text COLLATE pg_catalog."default",
    agency integer,
    zone_id text COLLATE pg_catalog."default",
    stop_url text COLLATE pg_catalog."default",
    location_type numeric,
    parent_station text COLLATE pg_catalog."default",
    stop_id text COLLATE pg_catalog."default"
TABLESPACE pg_default;

My current solution involves pandas and scipy kdtree and thought there must be a better way


First, you should add a geometry(point, 4326) field to your stops table, since you started with the other table as a geometry.

ALTER TABLE public.stops ADD COLUMN geom geometry(POINT,4326);
UPDATE public.stops SET geom =  ST_SetSRID(ST_Point(stop_lon, stop_lat),4326);

Make sure to have an index on the fields, for this case you better need a geography index:

CREATE INDEX stay_points_geog_idx ON public.stay_points USING GIST(cast(geom as geography));
CREATE INDEX stops_geog_idx ON public.stops USING GIST(cast(geom as geography));

Analyze your tables, to be sure that the planner knows about your indexes:

ANALYZE public.stops;
ANALYZE public.stay_points;

Now you can request the nearest-neighbor for your points, using postgis tool made for that (<->) and a CROSS JOIN LATERAL (don't forget to cast in geography):

                ST_Distance(geography(stay_points.geom), geography(stops.geom)) as distance
            FROM stops
            ORDER BY geography(stay_points.geom) <-> geography(stops.geom)
            LIMIT 1
        ) AS stops;

(untested code, but this should be the idea)

| improve this answer | |
  • Technically, you could do a LEFT JOIN LATERAL instead but the CROSS JOIN is common in this case because with that you can change the LIMIT 1 to a LIMIT10 for example to have for each stay_point the 10 nearest neighbors (and with a CROSS you will end up with 10 lines by stay_points in your result, whereas with the lateral still only the closest) – robin loche Oct 26 '18 at 9:58

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.