x cd /pub/data/sf12000/Tools;  %let pgm=aggsf1_moregns;  filename pgm "&pgm..sas";
/*-----------------------------------------------------------------------------------------------------------------
       This code aggregates sf12000 data to various Missouri regions, creating ph and pct data sets for:
         rpc  - Regional Planning Commissions
         ded  - Dept of Economic Development regions
         dot  - Dept of Transportaion districts
         uoe  - UM Outreach and Extension regions.

       Coded by John Blodgett,  OSEDA, U. of Missouri, blodgettj@umsystem.edu
       Under contract with the Missouri Census Data Center

      Revision History:
       3.22.02:  Coding begins.
 ------------------------------------------------------------------------------------------------------*/

 libname sf1  '/pub/data/sf12000';
 libname user '/tmp/scratch/user';    *<===Optional. Permissions on this directory limit who can use it;
 libname georef '/pub/data/georef';

 title 'sf12000 Aggregation Run: To Missouri Regions: rpcs, ded, dot and uoe';
 options msglevel=i;


 %macro atoI(list,x);
  %*--utility macro to take a list of variables containing a wild card character ("x", by default) and replacing it with the letters
      a, b, c, d, e, f, g, h, and I in order.  To generate the lists used in sf1 for the race/hispanic tables.--;
  %if &x eq %then %let x=x;   %*--default wild card character is "x"--;
  %sysfunc(translate(&list,a,&x))
  %sysfunc(translate(&list,b,&x))
  %sysfunc(translate(&list,c,&x))
  %sysfunc(translate(&list,d,&x))
  %sysfunc(translate(&list,e,&x))
  %sysfunc(translate(&list,f,&x))
  %sysfunc(translate(&list,g,&x))
  %sysfunc(translate(&list,h,&x))
  %sysfunc(translate(&list,I,&x))
  %mend atoI;

  /*<=========================Comment out=====================================
 data agginsf1(compress=yes);
  length sumlev $3 region $3;
  merge georef.mocogeos(keep=county ded_region rpc umx2 dot in=ingeos)
        sf1.moph(in=insf1 where=(sumlev='050'));
  by county;
  if not (ingeos and insf1) then do;
    put //'*******************PROBLEM!!!!!!!!!!!!!!!!!!!!!!: Data Missing on One of the files!!!!! '/
    county= ingeos= insf1=;
    abort abend;
    end;

 *---Output once for each region type--;
 keep sumlev region _numeric_;
    sumlev='ded';  region=ded_region;  output;
    sumlev='dot';  region=dot;   output;
    sumlev='umx';  region=umx2;  output;
    sumlev='rpc';  region=rpc;   output;
  run;

 proc sort data=agginsf1; by sumlev region; run;

 options compress=yes;    *<---to get the aggregated output set to be compressed--;

 %agg(aggin=agginsf1,aggout=aggoutsf1,
      aggby=sumlev region,
      aggvars=pop100 hu100 p1i1--h16Ii19 IntPtLon IntPtLat AreaSQMI,
      naggvars=,
      means=  IntPtLon IntPtLat p17i1  p33i1   p17ai1  p17bi1  p17ci1  p17di1  p17ei1  p17fi1 p17gi1  p17hi1  p17Ii1
                             p33ai1  p33bi1  p33ci1  p33di1  p33ei1  p33fi1 p33gi1  p33hi1  p33Ii1
                             h12i1   h12i2   h12i3
                             %atoI(h12xi1 h12xi2 h12xi3) ,
      mweights=p1i1     p1i1   p15i1  p31i1   p15ai1  p15bi1  p15ci1  p15di1  p15ei1  p15fi1 p15gi1  p15hi1  p15Ii1
                             p31ai1  p31bi1  p31ci1  p31di1  p31ei1  p31fi1 p31gi1  p31hi1  p31Ii1
                             h4i1    h4i2    h4i3
                             %atoI(h15xi1 h15xi2 h15xi10) ,
      debug=1,
      agglvls=1)
  run;
 ===========================End commented out================================ */

 data sf1.moregnsph(label='Missouri County-based Regions P & H Tables'  sortedby=SumLev Region) ;
  length sumlev $3 geocode region ded_Region $3 dot rpc umx2 $2  areaname $40 ;
  set aggoutsf1;
  geocode=region;
 select (sumlev);
   when('ded') do;
    ded_region=region;  areaname=put(ded_region,$dedrgn.);
    end;
  when('dot') do;
    dot=region;   areaname=put(dot,$modot.);
    end;
  when('rpc') do;
    rpc=region;  areaname=put(rpc,$rpcname.);
    end;
  when('umx') do;
    umx2=region;  areaname=umx2||' MU Extension Region';
    end;
  otherwise do;
    put sumlev= ' *****unexpectd sumlev code ----aborting step!!!!!!!!!!!!! ********************** ';
    abort abend;
    end;
  end;  *--select--*;

  drop _lvl_ _nag_ ;
 label ded_Region='Dept of Economic Development region' dot='MO Dept of Transp. District'
       rpc='Regional Planning Commission'  umx2='MU Extension Region';
 format intptlon 11.6 intptlat 10.6;

  *---Calculate all the age medians from the corresponding distributions.  Note that we do not have explicit distrib
      tables for total age, but have to create the totals by summing male and female counts in each case.--;
  array _p12t {23} _temporary_;
  array _p12w {23} _temporary_ (5, 5, 5, 3, 2,1,1,3,5,5,5,5,5,5,5,2,3,2,3,5,5,5,5);
  option  NOmprint;

  %macro domedian(r);
    %*--r is the race code for the table.  We do the same processing on corresponding p13 and p12 tables depending on
        the race code.  When r is blank we are working with the unqualified total-pop tables--;
    array _p12&r.m p12&r.i3-p12&r.i25;
    array _p12&r.f p12&r.i27-p12&r.i49;
    do _i_ = 1 to 23;
       _p12t{_i_}=_p12&r.m{_i_} + _p12&r.f{_i_};
       end;
    %*--Example: when r=a then the first macro invocation says to calculate the est median variable p13ai1 by using
        the data distribution values in array _p12t (data from table p12a, summed over sex) and using the array _p12w
        of constant "interval widths" data.   The second median invocation repeats this process but assigns result to
        p13ai2 and uses the _p12am data distribution array which is the set of male counts.  And the 3rd invocation
        creates p13ai3 as the female median by working with the array of females counts in _p12af.---;

    %median(mdn=p13&r.i1,intvals=_p12t,intsize=_p12w,samestep=1) ;    %*--medians for total pop-;
    %median(mdn=p13&r.i2,intvals=_p12&r.m,intsize=_p12w,samestep=1) ; %*--medians for males-;
    %median(mdn=p13&r.i3,intvals=_p12&r.f,intsize=_p12w,samestep=1) ; %*--medians for females-;
    %mend domedian;

  %domedian()   ;*--Assigns p13i1, p13i2 and p13i3  using data in distribution array p12-;
  %domedian(a)  ;*--Assigns p13ai1, p13ai2 and p13ai3  using data in distribution array p12a-;
  %domedian(b)  ;*--etc.--;
  %domedian(c)
  %domedian(d)
  %domedian(e)
  %domedian(f)
  %domedian(g)
  %domedian(h)
  %domedian(I)
  drop _haf -- _slot2;  *--temp vars used by median macro--;
   run;

 *---verify that the total pops agg to state total---;
proc means data=sf1.moregnsph(keep=sumlev pop100) sum n min max mean;
 class sumlev;
 var pop100;
 title 'Means of sf1.moregnsph, by SumLev';
 run;

*-----------------------Part 2: Aggregate the pct tables to the regions.-------------------------------------- ;

data agginpct;
  length SumLev Region $3;
  merge georef.mocogeos (in=inclist keep=county ded_region dot umx2 rpc)
        sf1.mopct(in=insf1 where=(sumlev='050')
            keep=sumlev county pop100 hu100 pct1i1 -- pct17Ii75);
  by county;
  if not (inclist and insf1) then  do;
    put //'*******Logic or Dataset Error: Not a 1-to-1 match between county level sf1 pct and mocogeos***** '/
    county= insf1=  inclist= ;
    abort abend;
    end;
  sumlev='ded';  region=ded_region;  output;
  sumlev='dot';  region=dot;   output;
  sumlev='umx';  region=umx2;  output;
  sumlev='rpc';  region=rpc;   output;
  keep sumlev region _numeric_;
  run;

proc sort data=agginpct;  by  sumlev region;  run;

*----Aggregate the pct table data by the districts----;


 %agg(aggin=agginpct,aggout=aggoutpct,
      aggby=sumlev region ,
      aggvars=pop100 hu100  pct1i1 -- pct17Ii75,
      naggvars=1920,
      dropvars=_lvl_ _nag_,
      agglvls=1)
  run;

 data sf1.moregnspct(label='Missouri County-based Regions PCT Tables' sortedby=SumLev Region);
  length sumlev $3 geocode region ded_Region $3 dot rpc umx2 $2  areaname $40 ;
  set aggoutpct;
  geocode=region;
  length sumlev $3 geocode region ded_region rpc  $3 dot umx2 $2 areaname $40 ;
  select (sumlev);
   when('ded') do;
    ded_region=region;  areaname=put(ded_region,$dedrgn.);
    end;
   when('dot') do;
     dot=region;   areaname=put(dot,$modot.);
     end;
   when('rpc') do;
     rpc=region;  areaname=put(rpc,$rpcname.);
     end;
   when('umx') do;
     umx2=region;  areaname=umx2||' MU Extension Region';
     end;
   otherwise do;
     put sumlev= ' *****unexpectd sumlev code ----aborting step!!!!!!!!!!!!! ********************** ';
     abort abend;
     end;
   end;  *--select--*;
 label ded_Region='Dept of Economic Development region' dot='MO Dept of Transp. District'
       rpc='Regional Planning Commission'  umx2='MU Extension Region';
  run;

proc means data=sf1.moregnspct(keep=sumlev pop100) sum n min max mean;
 class sumlev;
 var pop100;
 title2 'Means of sf1.moregnspct, by SumLev';
 run;

 %include sascode(notify);