From: Regina Obe Date: Sun, 17 Jul 2011 16:11:08 +0000 (+0000) Subject: enhancements to support reverse geocoding of highway locations. Also add reverse... X-Git-Tag: 2.0.0alpha1~1198 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=a1282f5e90784693a5e1a8fe18797ee340abcc6c;p=postgis enhancements to support reverse geocoding of highway locations. Also add reverse geocode regress tests git-svn-id: http://svn.osgeo.org/postgis/trunk@7647 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/extras/tiger_geocoder/tiger_2010/geocode/other_helper_functions.sql b/extras/tiger_geocoder/tiger_2010/geocode/other_helper_functions.sql index ebb3244a0..1bf5cfc27 100644 --- a/extras/tiger_geocoder/tiger_2010/geocode/other_helper_functions.sql +++ b/extras/tiger_geocoder/tiger_2010/geocode/other_helper_functions.sql @@ -116,6 +116,17 @@ SELECT 'CREATE INDEX idx_' || c.table_schema || '_' || c.table_name || '_l' || c FROM (SELECT table_name, table_schema FROM information_schema.tables WHERE table_type = 'BASE TABLE' AND table_name LIKE '%featnames' AND table_schema IN('tiger','tiger_data')) As t INNER JOIN (SELECT * FROM information_schema.columns WHERE column_name IN('fullname') ) AS c + ON (t.table_name = c.table_name AND t.table_schema = c.table_schema) + LEFT JOIN pg_catalog.pg_indexes i ON + (i.tablename = c.table_name AND i.schemaname = c.table_schema + AND indexdef LIKE '%btree%(%lower%' || c.column_name || ')%varchar_pattern_ops%') +WHERE i.tablename IS NULL +-- var_ops mtfcc -- +UNION ALL +SELECT 'CREATE INDEX idx_' || c.table_schema || '_' || c.table_name || '_' || c.column_name || '_var_ops' || ' ON ' || c.table_schema || '.' || c.table_name || ' USING btree(' || c.column_name || ' varchar_pattern_ops);' As index +FROM (SELECT table_name, table_schema FROM + information_schema.tables WHERE table_type = 'BASE TABLE' AND table_name LIKE '%featnames' or table_name LIKE '%edges' AND table_schema IN('tiger','tiger_data')) As t INNER JOIN + (SELECT * FROM information_schema.columns WHERE column_name IN('mtfcc') ) AS c ON (t.table_name = c.table_name AND t.table_schema = c.table_schema) LEFT JOIN pg_catalog.pg_indexes i ON (i.tablename = c.table_name AND i.schemaname = c.table_schema diff --git a/extras/tiger_geocoder/tiger_2010/geocode/reverse_geocode.sql b/extras/tiger_geocoder/tiger_2010/geocode/reverse_geocode.sql index 95463c292..f30f08dab 100644 --- a/extras/tiger_geocoder/tiger_2010/geocode/reverse_geocode.sql +++ b/extras/tiger_geocoder/tiger_2010/geocode/reverse_geocode.sql @@ -14,15 +14,19 @@ $$ DECLARE var_redge RECORD; var_state text := NULL; + var_stusps text := NULL; var_countyfp text := NULL; var_addy NORM_ADDY; - var_strnum varchar; + var_addy_alt NORM_ADDY; var_nstrnum numeric(10); var_primary_line geometry := NULL; var_primary_dist numeric(10,2) ; var_pt geometry; + var_place varchar; var_stmt text; var_debug boolean = false; + var_zip varchar := NULL; + var_primary_fullname varchar := ''; BEGIN --$Id$ IF pt IS NULL THEN @@ -42,7 +46,7 @@ BEGIN IF var_debug THEN RAISE NOTICE 'Get matching states start: %', clock_timestamp(); END IF; - SELECT statefp INTO var_state FROM state WHERE ST_Intersects(the_geom, var_pt) LIMIT 1; + SELECT statefp, stusps INTO var_state, var_stusps FROM state WHERE ST_Intersects(the_geom, var_pt) LIMIT 1; IF var_debug THEN RAISE NOTICE 'Get matching states end: % - %', var_state, clock_timestamp(); END IF; @@ -53,7 +57,18 @@ BEGIN IF var_debug THEN RAISE NOTICE 'Get matching counties start: %', clock_timestamp(); END IF; + -- locate county SELECT countyfp INTO var_countyfp FROM county WHERE statefp = var_state AND ST_Intersects(the_geom, var_pt) LIMIT 1; + + --locate zip + SELECT zcta5ce INTO var_zip FROM zcta5 WHERE statefp = var_state AND ST_Intersects(the_geom, var_pt) LIMIT 1; + + -- lcoate city + IF var_zip > '' THEN + var_addy.zip := var_zip ; + SELECT z.place INTO var_place FROM zip_state_loc AS z WHERE z.zip = var_zip and z.statefp = var_state LIMIT 1; + var_addy.location := var_place; + END IF; IF var_debug THEN RAISE NOTICE 'Get matching counties end: % - %',var_countyfp, clock_timestamp(); @@ -63,6 +78,8 @@ BEGIN RETURN; END IF; + var_addy.stateAbbrev = var_stusps; + -- Find the street edges that this point is closest to with tolerance of 0.005 but only consider the edge if the point is contained in the right or left face -- Then order addresses by proximity to road IF var_debug THEN @@ -70,37 +87,41 @@ BEGIN END IF; var_stmt := ' - WITH + WITH ref AS ( + SELECT ' || quote_literal(var_pt::text) || '::geometry As ref_geom ) , f AS - ( SELECT * FROM faces + ( SELECT * FROM faces CROSS JOIN ref WHERE statefp = ' || quote_literal(var_state) || ' AND countyfp = ' || quote_literal(var_countyfp) || ' - AND ST_DWithin(faces.the_geom, ' || quote_literal(var_pt::text) || '::geometry,0.0001) ), + AND ST_Intersects(faces.the_geom, ref_geom) + ), e AS - ( SELECT DISTINCT ON (edges.tlid) edges.* - FROM edges INNER JOIN f ON (edges.tfidr = f.tfid OR edges.tfidl = f.tfid) + ( SELECT edges.* , CASE WHEN edges.tfidr = f.tfid THEN ''R'' WHEN edges.tfidl = f.tfid THEN ''L'' ELSE NULL END::varchar As eside + FROM edges INNER JOIN f ON (f.statefp = edges.statefp AND (edges.tfidr = f.tfid OR edges.tfidl = f.tfid)) + CROSS JOIN ref WHERE edges.statefp = ' || quote_literal(var_state) || ' AND edges.countyfp = ' || quote_literal(var_countyfp) || ' - AND ST_DWithin(edges.the_geom, ' || quote_literal(var_pt::text) || '::geometry, 0.05) AND (edges.mtfcc BETWEEN ''R'' AND ''T'') --only consider streets and roads + AND ST_DWithin(edges.the_geom, ref.ref_geom, 0.01) AND (edges.mtfcc LIKE ''S%'') --only consider streets and roads ), - a AS (SELECT addr.* FROM addr INNER JOIN e ON addr.tlid = e.tlid - WHERE addr.statefp = ' || quote_literal(var_state) || ') + a AS (SELECT addr.* FROM addr INNER JOIN e ON (addr.statefp = e.statefp AND addr.tlid = e.tlid) + WHERE addr.statefp = ' || quote_literal(var_state) || '), + n AS (SELECT featnames.* FROM featnames INNER JOIN e ON(featnames.statefp = e.statefp AND featnames.tlid = e.tlid) + WHERE featnames.statefp = ' || quote_literal(var_state) || ' AND e.mtfcc LIKE ''S%'' ) SELECT * - FROM (SELECT DISTINCT ON(fullname) foo.fullname, foo.stusps, foo.zip, - (SELECT z.place FROM zip_state_loc AS z WHERE z.zip = foo.zip and z.statefp = ' || quote_literal(var_state) || ' LIMIT 1) As place, foo.center_pt, - side, to_number(fromhn, ''999999'') As fromhn, to_number(tohn, ''999999'') As tohn, ST_GeometryN(ST_Multi(line),1) As line, foo.dist + FROM (SELECT DISTINCT ON(tlid,eside) foo.fullname, foo.streetname, foo.streettypeabbrev, foo.zip, foo.center_pt, + eside, to_number(fromhn, ''999999'') As fromhn, to_number(tohn, ''999999'') As tohn, ST_GeometryN(ST_Multi(line),1) As line, + ST_Distance_Sphere(foo.line, ' || quote_literal(var_pt::text) || '::geometry) As dist FROM - (SELECT e.the_geom As line, e.fullname, a.zip, s.abbrev As stusps, ST_ClosestPoint(e.the_geom,' || quote_literal(var_pt::text) || '::geometry) As center_pt, e.statefp, a.side, a.fromhn, a.tohn, - ST_Distance_Sphere(e.the_geom, ' || quote_literal(var_pt::text) || '::geometry) As dist + (SELECT e.tlid, e.the_geom As line, COALESCE(n.fullname,e.fullname) As fullname, COALESCE(n.prequalabr || '' '','''') || n.name AS streetname, n.predirabrv, COALESCE(suftypabrv, pretypabrv) As streettypeabbrev, + n.sufdirabrv, a.zip, ST_ClosestPoint(e.the_geom,' || quote_literal(var_pt::text) || '::geometry) As center_pt, e.eside, a.fromhn, a.tohn , + ST_Distance_Sphere(e.the_geom, ' || quote_literal(var_pt::text) || '::geometry) As dist FROM e - INNER JOIN (SELECT * FROM state_lookup WHERE statefp = ' || quote_literal(var_state) || ' ) As s ON (e.statefp = s.statefp ) - LEFT JOIN f As fl ON (e.tfidl = fl.tfid) - LEFT JOIN f As fr ON (e.tfidr = fr.tfid) LEFT JOIN a ON - ( a.tlid = e.tlid AND ( ST_Covers(fl.the_geom, ' || quote_literal(var_pt::text) || '::geometry) AND a.side = ''L'') OR ( ST_Covers(fr.the_geom, ' || quote_literal(var_pt::text) || '::geometry) AND a.side = ''R'' ) ) - -- INNER JOIN zip_state_loc As z ON (a.statefp = z.statefp AND a.zip = z.zip) /** really slow with this join **/ - ORDER BY ST_Distance(e.the_geom, ' || quote_literal(var_pt::text) || '::geometry) LIMIT 4) As foo - WHERE dist < 200 --less than 150 m - ORDER BY foo.fullname, foo.dist) As f ORDER BY f.dist '; + ( a.tlid = e.tlid + AND a.side = e.eside + ) + LEFT JOIN n ON (n.statefp = e.statefp AND n.tlid = e.tlid) + ORDER BY dist LIMIT 50 ) As foo + ORDER BY foo.tlid, foo.eside, foo.fullname ASC NULLS LAST, dist LIMIT 15) As f ORDER BY f.dist '; IF var_debug = true THEN RAISE NOTICE 'Statement 1: %', var_stmt; @@ -122,46 +143,80 @@ BEGIN ORDER BY ST_Distance(e.the_geom, var_pt) LIMIT 4) As foo WHERE dist < 150 --less than 150 m ORDER BY foo.fullname, foo.dist) As f ORDER BY f.dist LOOP **/ - FOR var_redge IN EXECUTE var_stmt LOOP - IF var_debug THEN - RAISE NOTICE 'Get matching edges loop: %,%', var_primary_line, clock_timestamp(); - END IF; - IF var_primary_line IS NULL THEN --this is the first time in the loop and our primary guess - var_primary_line := var_redge.line; - var_primary_dist := var_redge.dist; - END IF; - -- We only consider other edges as matches if they intersect our primary edge -- that would mean we are at a corner place - IF ST_Intersects(var_redge.line, var_primary_line) THEN - intpt := array_append(intpt,var_redge.center_pt); - IF var_redge.fullname IS NOT NULL or var_redge.place IS NOT NULL THEN - street := array_append(street, (CASE WHEN include_strnum_range THEN COALESCE(var_redge.fromhn::varchar, '')::varchar || ' - ' || COALESCE(var_redge.tohn::varchar,'')::varchar || ' '::varchar ELSE '' END::varchar || COALESCE(var_redge.fullname::varchar,''))::varchar); - --interploate the number -- note that if fromhn > tohn we will be subtracting which is what we want - -- We only consider differential distances are reeally close from our primary pt - IF var_redge.dist < var_primary_dist*1.1 THEN - var_nstrnum := (var_redge.fromhn + ST_Line_Locate_Point(var_redge.line, var_pt)*(var_redge.tohn - var_redge.fromhn))::numeric(10); - -- The odd even street number side of street rule - IF (var_nstrnum % 2) != (var_redge.tohn % 2) THEN - var_nstrnum := CASE WHEN var_nstrnum + 1 NOT BETWEEN var_redge.fromhn AND var_redge.tohn THEN var_nstrnum - 1 ELSE var_nstrnum + 1 END; - END IF; - var_strnum := var_nstrnum::varchar; - IF var_redge.fullname IS NULL THEN --not a full address so can't be fully normalized - var_addy := NULL; - var_addy.address = var_strnum; - var_addy.location = var_redge.place; - var_addy.stateAbbrev = var_redge.stusps; - var_addy.zip = var_redge.zip; - ELSE - var_addy := normalize_address( COALESCE(var_strnum::varchar || ' ', '') || COALESCE(var_redge.fullname || ', ') || var_redge.place || ', ' || var_redge.stusps || ' ' || COALESCE(var_redge.zip,'')); - END IF; - addy := array_append(addy, var_addy); - END IF; - END IF; - END IF; - END LOOP; + + FOR var_redge IN EXECUTE var_stmt LOOP + IF var_debug THEN + RAISE NOTICE 'Start Get matching edges loop: %,%', var_primary_line, clock_timestamp(); + END IF; + IF var_primary_line IS NULL THEN --this is the first time in the loop and our primary guess + var_primary_line := var_redge.line; + var_primary_dist := var_redge.dist; + END IF; - IF var_debug THEN - RAISE NOTICE 'Get matching edges loop: %', clock_timestamp(); + IF var_redge.fullname IS NOT NULL AND COALESCE(var_primary_fullname,'') = '' THEN -- this is the first non-blank name we are hitting grab info + var_primary_fullname := var_redge.fullname; + var_addy.streetname = var_redge.streetname; + var_addy.streettypeabbrev := var_redge.streettypeabbrev; + END IF; + + IF ST_Intersects(var_redge.line, var_primary_line) THEN + var_addy.streetname := var_redge.streetname; + + var_addy.streettypeabbrev := var_redge.streettypeabbrev; + var_addy.address := var_nstrnum; + IF var_redge.fromhn IS NOT NULL THEN + --interpolate the number -- note that if fromhn > tohn we will be subtracting which is what we want + var_nstrnum := (var_redge.fromhn + ST_Line_Locate_Point(var_redge.line, var_pt)*(var_redge.tohn - var_redge.fromhn))::numeric(10); + -- The odd even street number side of street rule + IF (var_nstrnum % 2) != (var_redge.tohn % 2) THEN + var_nstrnum := CASE WHEN var_nstrnum + 1 NOT BETWEEN var_redge.fromhn AND var_redge.tohn THEN var_nstrnum - 1 ELSE var_nstrnum + 1 END; + END IF; + var_addy.address := var_nstrnum; + END IF; + IF var_redge.zip > '' AND COALESCE(var_addy.zip,'') <> var_redge.zip THEN + var_addy.zip := var_redge.zip; + SELECT z.place INTO var_place FROM zip_state_loc AS z WHERE z.zip = var_redge.zip and z.statefp = var_state LIMIT 1; + var_addy.location := var_place; + END IF; + + -- This is a cross streets - only add if not the primary adress street + IF var_redge.fullname > '' AND var_redge.fullname <> var_primary_fullname THEN + street := array_append(street, (CASE WHEN include_strnum_range THEN COALESCE(var_redge.fromhn::varchar, '')::varchar || COALESCE(' - ' || var_redge.tohn::varchar,'')::varchar || ' '::varchar ELSE '' END::varchar || COALESCE(var_redge.fullname::varchar,''))::varchar); + END IF; + + -- consider this a potential address + IF (var_redge.dist < var_primary_dist*1.1 OR var_redge.dist < 20) THEN + -- We only consider this a possible address if it is really close to our point + + intpt := array_append(intpt,var_redge.center_pt); + addy := array_append(addy, var_addy); + -- note that ramps don't have names or addresses but they connect at the edge of a range + -- so for ramps the address of connecting is still useful + + -- Determine city if zip is different from previous + + END IF; + END IF; + END LOOP; + + -- not matching roads or streets, just return basic info + IF NOT FOUND THEN + addy := array_append(addy,var_addy); + END IF; + + IF var_addy.streetname > '' AND addy[1].streetname IS NULL THEN + --there were no intersecting addresses with names, just pick the best candidate - the match is proably an offshoot of some sort + var_addy_alt := addy[1]; + var_addy_alt.streetname := var_addy.streetname; + var_addy_alt.streettypeabbrev := var_addy.streettypeabbrev; + addy[1] := var_addy_alt; + END IF; + IF var_debug THEN + RAISE NOTICE 'End Get matching edges loop: %', clock_timestamp(); END IF; + + + RETURN; END; diff --git a/extras/tiger_geocoder/tiger_2010/regress/regress.sql b/extras/tiger_geocoder/tiger_2010/regress/regress.sql index 5a2b7095d..3b37a907f 100644 --- a/extras/tiger_geocoder/tiger_2010/regress/regress.sql +++ b/extras/tiger_geocoder/tiger_2010/regress/regress.sql @@ -2,4 +2,6 @@ \o normalize_address_regress.out \i normalize_address_regress.sql \o geocode_regress.out -\i geocode_regress.sql \ No newline at end of file +\i geocode_regress.sql +\o reverse_geocode_regress.out +\i reverse_geocode_regress.sql \ No newline at end of file diff --git a/extras/tiger_geocoder/tiger_2010/regress/reverse_geocode_regress b/extras/tiger_geocoder/tiger_2010/regress/reverse_geocode_regress new file mode 100644 index 000000000..9a4b62335 --- /dev/null +++ b/extras/tiger_geocoder/tiger_2010/regress/reverse_geocode_regress @@ -0,0 +1,4 @@ +I- 90, Wellesley, MA 02493|{"(,,90,I-,,,Wellesley,MA,02493,)"} +I- 90, Worcester, MA 01501|{"(,,90,I-,,,Worcester,MA,01501,)"} +158 Washington St, Boston, MA 02108|{"(158,,Washington,St,,,Boston,MA,02108,)"} +32 Capen St, Medford, MA 02155|{"(32,,Capen,St,,,Medford,MA,02155,)","(3,,Edison,Ave,,,Medford,MA,02155,)"} diff --git a/extras/tiger_geocoder/tiger_2010/regress/reverse_geocode_regress.sql b/extras/tiger_geocoder/tiger_2010/regress/reverse_geocode_regress.sql new file mode 100644 index 000000000..50ab3b7b9 --- /dev/null +++ b/extras/tiger_geocoder/tiger_2010/regress/reverse_geocode_regress.sql @@ -0,0 +1,6 @@ +\timing +SELECT pprint_addy(addy[1]), addy FROM reverse_geocode(ST_Point(-71.27593,42.33891)); -- I 90 Exit 14 +SELECT pprint_addy(addy[1]), addy FROM reverse_geocode(ST_Point(-71.85335,42.19262)); -- I 90 Exit 10, Worcester MA +SELECT pprint_addy(addy[1]), addy FROM reverse_geocode(ST_Point(-71.057811,42.358274)); -- 1 Devonshire Place (washington st area) +SELECT pprint_addy(addy[1]), addy FROM reverse_geocode(ST_Point(-71.123848,42.41115)); --30 capen +\timing \ No newline at end of file