/*-----
 * group: Miscellaneous
 * purpose: Return list of zipcodes within X miles of another
 * notes: <A HREF="50.htm">The $50 dollar macro</A> !
 */

%macro getzips (data=tables.zipcode, findzip=?, xmiles=, out=work.gotzips);

%local data findzip xmiles out;

%* Richard A. DeVenezia, Copyright 2000
%* 9949 East Steuben Road
%* Remsen, NY 13438
%* contact: http://www.devenezia.com/contact.php
%*
%* Coded:  9/25/00
%* Update: 4/ 2/04 fixed problem (reported by Zhonghe Li) with %if after first proc sql
%*
%* Purpose:
%*   Determine zipcodes within (<=) xmiles of zipcode findzip
%*   Distance is computed as 'spherical distance' using the Haversine Formula
%*
%* Related Topic:
%*   Calculate distance between two (2) points given latitude and longitude
%*
%* Arguments:
%*   data    - table of zip codes and 'centroid' location in latitude and longitude
%*   findzip - zipcode to which distances are computed
%*   xmiles  - zipcodes within this radius (miles) of zipcode are to be returned
%*   out     - output table that will contain only zipcodes
%*
%* Requirements:
%* zipdata
%*   zipcode column (char) must be named ZIPCODE
%*   latitude column (number) must be named LATITUDE and units degrees
%*   longitude column (number) must be named LONGTUDE and units degrees
%*
%* Caveats:
%*   No assumptions about centroid locations being projected or unprojected.
%*   Distance between points is 'surface distance', see
%*   http://earth.uni-muenster.de/~eicksch/GMT-Help/msg00147.html
%*   or search the web for "surface distance between points on map"
%*   or consult the GIS FAQ (http://www.census.gov/geo/www/faq-index.html)
%*   If a validation error occurs, the output data set is not touched
%*
%* Formula:
%*   Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
%*   Sky and Telescope, vol. 68, no. 2, 1984, p. 159):
%*   dlon = lon2 - lon1
%*   dlat = lat2 - lat1
%*   a = sin^2(dlat/2) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
%*   c = 2 * arcsin(min(1,sqrt(a)))
%*   d = R * c
%*   where
%*   lat and lon are latitude and longitude in radian respectively
%*   radians = degrees * pi / 180, (pi/180) = 0.017453293
%*   R = 3956, radius of Earth in miles
%*;

  %local radius d2r;

  %let radius = 3956;
  %let d2r = 0.017453293; %* radians / degree ;

  %* validate zipcode data set;

  %local zipCol latCol longCol;

  %let zipCol = ZIPCODE;
  %let latCol = LATITUDE;
  %let longCol= LONGTUDE;

  %local error dsid vn1 vn2 vn3 vt1 vt2 vt3;

  %let error = 0;

  %let dsid = %sysfunc (open (&data));
  %if &dsid > 0 %then %do;
    %let vn1 = %sysfunc (varnum (&dsid, &zipCol));
    %let vn2 = %sysfunc (varnum (&dsid, &latCol));
    %let vn3 = %sysfunc (varnum (&dsid, &longCol));

    %if &vn1 = 0 %then %do;
      %put ERROR: &data does not contain column &zipCol..;
      %let error = 1;
    %end;
    %if &vn2 = 0 %then %do;
      %put ERROR: &data does not contain column &latCol..;
      %let error = 1;
    %end;
    %if &vn3 = 0 %then %do;
      %put ERROR: &data does not contain column &longCol..;
      %let error = 1;
    %end;
    %if &error %then %goto ByeBye;

    %let vt1 = %sysfunc (vartype (&dsid, &vn1));
    %let vt2 = %sysfunc (vartype (&dsid, &vn2));
    %let vt3 = %sysfunc (vartype (&dsid, &vn3));

    %if &vt1 ne C %then %do;
      %put ERROR: column &zipCol is of type &vt1 and should be character.;
      %let error = 1;
    %end;
    %if &vt2 ne N %then %do;
      %put ERROR: column &latCol is of type &vt2 and should be numeric.;
      %let error = 1;
    %end;
    %if &vt3 ne N %then %do;
      %put ERROR: column &longCol is of type &vt3 and should be numeric.;
      %let error = 1;
    %end;
    %if &error %then %goto ByeBye;

    %let dsid = %sysfunc (close (&dsid));
  %end;
  %else %do;
    %put ERROR: &data could not be opened.;
    %let error = 1;
  %end;
  %if &error %then %goto ByeBye;

  %* validate findzip;

  %local zip1;

  %let zip1 = %sysfunc (dequote (&findzip));

  %if %length (&zip1) = %length (&findzip) %then %do;
    %* findzip was not dequoted, assume only a number was passed;
    %if &zip1 ne ? %then %do;
      %let zip1 = %sysfunc (putn (&zip1, z5.));
      %let zip1 = "&zip1";
      %put WARNING: Argument FINDZIP=&findzip was converted to zipcode &zip1..;
    %end;
    %else %do;
      %put ERROR: Argument FINDZIP= was not specified.;
      %goto ByeBye;
    %end;
  %end;
  %else %do;
    %* findzip was dequoted, use original argument;
    %let zip1 = &findzip;
  %end;

  %* obtain the latitude and longitude of findzip (and convert to radians);

  %local lat1r long1r;

  %let lat1r = -1;
  %let long1r = -1;

  Proc SQL noprint;
    select min (&latCol) * &d2r, min (&longCol) * &d2r
    into :lat1r, :long1r
    from &data
    where &zipCol = &findzip
    having count(*) > 0;
  quit;

  %if %sysevalf(&lat1r = -1) or %sysevalf(&long1r = -1) %then %do;
    %put ERROR: &zipCol &findzip was not found in &data;
    %goto ByeBye;
  %end;

  %* select the zipcodes whose sphere surface distance is <= from zip1;

  Proc SQL;
    create table &out (keep=&zipCol) as
    select
      &zipCol
    , &latCol  * &d2r as lat2r
    , &longCol * &d2r - &long1r as dLon
    , calculated lat2r - &lat1r  as dLat
    , sin(calculated dLat/2)**2 as a1
    , cos(&lat1r) as a2
    , cos(calculated lat2r) as a3
    , sin(calculated dLon/2)**2 as a4
    ,   calculated a1
      + calculated a2 * calculated a3 * calculated a4 as a
    , 2 * arsin (min (1, sqrt(calculated a))) as c
    , &radius * calculated c as distance
    from &data
    where
      calculated distance <= &xmiles
    ;
    ;
  quit;

%ByeBye:
  %* make sure zipcode data set is closed when error exiting;
  %if &dsid %then
    %let dsid = %sysfunc (close (&dsid));

%mend getzips;

/*
%getzips (data=ctzip, findzip = "06013", xmiles = 25 );
*/