First draft of reverse geocde currently just gets the street need to interpolate...

git-svn-id: http://svn.osgeo.org/postgis/trunk@6817 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/extras/tiger_geocoder/tiger_2010/geocode/reverse_geocode.sql b/extras/tiger_geocoder/tiger_2010/geocode/reverse_geocode.sql
new file mode 100644 (file)
index 0000000..62a7830
--- /dev/null
+++ b/extras/tiger_geocoder/tiger_2010/geocode/reverse_geocode.sql
@@ -0,0 +1,65 @@
+--$Id$\r
+ /*** \r
+ * \r
+ * Copyright (C) 2011 Regina Obe and Leo Hsu (Paragon Corporation)\r
+ **/\r
+-- This function given a point try to determine the approximate street address (norm_addy form)\r
+-- and array of cross streets, as well as interpolated points along the streets\r
+-- Use case example: SELECT pprint_addy(r.addy[1]) As st1, array_to_string(r.street, ',') FROM reverse_geocode(ST_GeomFromText('POINT(-71.057811 42.358274)',4269)) As r;\r
+set search_path=tiger,public;\r
+CREATE OR REPLACE FUNCTION reverse_geocode(\r
+    IN pt geometry,\r
+    IN include_strnum_range boolean,\r
+    OUT intpt geometry[],\r
+    OUT addy NORM_ADDY[],\r
+    OUT street text[]\r
+) RETURNS RECORD\r
+AS $_$\r
+DECLARE\r
+  var_redge RECORD;\r
+  var_states text[];\r
+  var_addy NORM_ADDY;\r
+BEGIN\r
+  IF pt IS NULL THEN\r
+    RETURN;\r
+  END IF;\r
+  -- Determine state tables to check \r
+  -- this is needed to take advantage of constraint exclusion\r
+  var_states := ARRAY(SELECT statefp FROM state WHERE ST_Intersects(the_geom, pt) );\r
+  IF array_upper(var_states, 1) IS NULL THEN\r
+  -- We don't have any data for this state\r
+       RETURN;\r
+  END IF;\r
+  \r
+  -- Find the street edges that this point is closest to with tolerance of 500 meters (width of road)\r
+  FOR var_redge IN\r
+        SELECT DISTINCT ON(fullname) foo.gid, foo.fullname, foo.stusps, foo.zip, (SELECT z.place FROM zip_state_loc AS z WHERE z.zip = foo.zip and z.statefp = foo.statefp) As place, foo.intpt \r
+        FROM \r
+          (SELECT e.gid, e.fullname, a.zip, s.stusps, ST_ClosestPoint(e.the_geom, pt) As intpt, e.statefp\r
+                       FROM edges AS e INNER JOIN state As s ON (e.statefp = s.statefp AND s.statefp = ANY(var_states) )\r
+                               INNER JOIN addr As a ON (e.tlid = a.tlid AND e.statefp = a.statefp) -- INNER JOIN zip_state_loc As z ON (a.statefp =  z.statefp AND a.zip = z.zip) /** really slow with this join **/\r
+                       WHERE e.statefp = ANY(var_states) AND a.statefp = ANY(var_states) AND ST_DWithin(e.the_geom, pt, 0.005)\r
+                       ORDER BY ST_Distance(e.the_geom, pt) LIMIT 4) As foo \r
+                       WHERE ST_Distance_Sphere(foo.intpt, pt) < 150 --less than 150 m\r
+                       ORDER BY foo.fullname, ST_Distance_Sphere(foo.intpt, pt) LOOP\r
+       intpt := array_append(intpt,var_redge.intpt); \r
+       IF var_redge.fullname IS NOT NULL THEN\r
+                       street := array_append(street,var_redge.fullname::text);\r
+                       var_addy := normalize_address(var_redge.fullname || ', ' || var_redge.place || ', ' || var_redge.stusps || ' ' || var_redge.zip);\r
+                       addy := array_append(addy, var_addy);\r
+       END IF;\r
+  END LOOP;\r
+                       \r
+  RETURN;   \r
+END;\r
+$_$ LANGUAGE plpgsql;\r
+\r
+CREATE OR REPLACE FUNCTION reverse_geocode(IN pt geometry, OUT intpt geometry[],\r
+    OUT addy NORM_ADDY[],\r
+    OUT street text[]) RETURNS RECORD\r
+AS \r
+$$\r
+-- default to not including street range in cross streets\r
+SELECT reverse_geocode($1,false);\r
+$$\r
+language sql;
\ No newline at end of file