[ create a new paste ] login | about

Link: http://codepad.org/jt9WVoi7    [ raw code | fork ]

Plain Text, pasted on Oct 24:
CREATE OR REPLACE FUNCTION imt_rgeo_countyzip(geometry)
  RETURNS rgeo_out AS
$BODY$
/*
 * ROW(county, state, fips, zipcode, ct_dist_m, zp_dist_m) = rgeo(point);
 * ROW(county, state, fips, zipcode, ct_dist_m, zp_dist_m) = rgeo(long, lat);
 * 
 * This function performs reverse geocoding on county and zcta5 tables. 
 * It returns a rgeo_result record for the closest object, or
 * NULL if it fails.
 *
 * We are ignoring the fact the degrees longitude are shortened as we
 * move toword the poles.
 *
 * Algorithm:
 * We create a small radius (0.0013 degrees, .0897 miles) and search
 * for the closest object within that radius and compute the results.
 * If we fail to find any results in that radius we double it and search
 * again until we exceed a one degree radius which is roughtly 69*69 sq-miles.
 *
 */

DECLARE
    pnt ALIAS FOR $1;
    radius FLOAT;
    rr RECORD;
    zip RECORD;
    ret rgeo_out;
    max_radius FLOAT := 0.5;

BEGIN
    IF GeometryType(pnt) != 'POINT' THEN
        RETURN NULL;
    END IF;
raise notice 'calling imt_rgeo_countyzip(%,%)', st_x(pnt), st_y(pnt);

    radius := 0.0013;
    LOOP
        SELECT cname
             , stusps
             , stccc
             , ST_Distance_Spheroid(pnt, the_geom, 'SPHEROID["WGS 84",6378137,298.257223563]') as dist
            INTO rr
            FROM data.county
            WHERE st_dwithin(pnt, the_geom, radius)
            ORDER BY distance(pnt, the_geom) ASC limit 1;
        IF FOUND THEN
raise notice 'rr: %, radius: %', rr, radius;
            EXIT;
        END IF;
        radius := radius * 2.0;
        IF radius > max_radius THEN
	    rr.cname  = null;
	    rr.stusps = null;
	    rr.stccc  = null;
	    rr.dist   = null;
            EXIT;
        END IF;
    END LOOP;

    -- get the zipcode ignoring nnnHH and nnnXX
    radius := 0.0013;
    LOOP
        SELECT zipcode
             , ST_Distance_Spheroid(pnt, the_geom, 'SPHEROID["WGS 84",6378137,298.257223563]') as dist 
             INTO zip
             FROM data.zcta5
            WHERE st_dwithin(pnt, the_geom, radius)
              and zipcode not like '___HH' and zipcode not like '___XX'
            ORDER BY distance(pnt, the_geom), zipcode ASC limit 1;
        IF FOUND THEN
raise notice 'zip: %, radius: %', zip, radius;
            ret := ROW(rr.cname, rr.stusps, rr.stccc, zip.zipcode, rr.dist, zip.dist);
            RETURN ret;
        END IF;
        radius := radius * 2.0;
        IF radius > max_radius THEN
                ret := ROW(rr.cname, rr.stusps, rr.stccc, null, rr.dist, null);
                RETURN ret;
        END IF;
    END LOOP;
END
$BODY$
  LANGUAGE plpgsql STABLE
  COST 100;



Create a new paste based on this one


Comments: