/*--------------------------------------------------------------------*
 | Purpose:     Use PROC SGMAP to display wind power generation 
 |              capacities in the U.S. 
 |
 | Version:     SAS 9.40M5 or later
 |
 | Data Source: U.S. Energy Information Administration (Form EIA-860)
 |
 | Usage:       1) Download 2016 ZIP file from:
 |                    https://www.eia.gov/electricity/data/eia860/
 |              2) Extract these files to your local system:
 |                    2___Plant_Y2016.xlsx
 |                    3_2_Wind_Y2016.xlsx
 |              3) Set PATH variable to that extracted file location.
 |              4) Submit program to import files and create maps.
 |
 | Assistance:  communities.sas.com
 *--------------------------------------------------------------------*/

/*--- Set PATH variable to location of unzipped files from EIA, */
%let path=/folders/myfolders;

/*--- Import Excel file containing wind turbine parameters at all sites. */
proc import file="&path./3_2_Wind_Y2016.xlsx"
            out=turbine_parameters
            dbms=xlsx replace;
  /*--- Use Excel column headings as SAS variable names. */ 
  getnames=yes;
  /*--- Specify Excel sheet name and data range to import. */
  range='Operable$A2:AP1200';
run;

/*--- Sort the turbine data so we can sum the generation capacity
      at each site and later combine it with the plant data.
      We do not need all the variables from the imported files.
      Use 'keep' option to select the needed variables. */
proc sort data=turbine_parameters(keep=Utility_ID Plant_Code 
                                       Nameplate_Capacity__MW_
                                       Number_of_Turbines 
                                       Turbine_Hub_Height__Feet_);
  /*--- Specify variables to sort data by. We need to group the 
        turbine observations by utility plant owner and plant ID. */
  by Utility_ID Plant_Code;
run;

/*--- Data can contain multiple turbines at same site.
      Sum the turbine capacities for each site.  */
proc summary data=turbine_parameters noprint;
  /*--- Specify the values to sum. */
  var Nameplate_Capacity__MW_;
  /*--- Sum the capacities for each plant owner and each site. */
  by Utility_ID Plant_Code;
  /*--- Specify data set for the summary output, drop variables not
        needed and specify name of the variable (MW_total) to contain
        the summed capacities. */
  output out=turbines_by_site(drop=_type_ _freq_) sum=MW_total;
run;

/*--- Turbine data does not include locations. Import file that
      contains latitude and longitude for each plant. */
proc import file="&path./2___Plant_Y2016.xlsx"
            out=power_plants
            dbms=xlsx replace;
  getnames=yes;
  range='Plant$A2:AP10000';
run;

/*--- Sort plant locations to combine with turbine data.
      Note that the 'by' statement lists the same variables used when
      the turbine data was sorted above. */
proc sort data=power_plants
          (keep=Utility_ID Utility_Name Plant_Code Latitude Longitude State);
  by Utility_ID Plant_Code;
run;

/*--- Combine power plant locations with turbine data. Write observations into two
      data sets. First contains all wind power sites in the U.S. with more than 
      200 MW of capacity. Second contains all sites in the New England states. */
data US_GT_200MW(label='U.S. Energy Information Administration, 2016 Wind Power Generation >200 MW')
     New_England_All(label='U.S. Energy Information Administration, 2016 Wind Power Generation in New England');
  /*--- Combine the turbine data with the plant location data. 
        Tag turbine data with 'a' so we can use it as a filter below. */
  merge turbines_by_site (in=a) power_plants;
  /*--- Specify values to use when combining the observations. 
        We want to combine values which have the same utility owner 
        and at the same power plant. */
  by Utility_ID Plant_Code;
  /*--- Add a label to explicitly describe the 'MW_total' values. */
  label MW_total='Total turbine capacity at site (MW)';
  /*--- SGMAP requires cooridinates be numeric variables but the import
        process made them character variables. Convert them to numerics. */
  lat_num  = input( latitude, best. );
  long_num = input( longitude, best. );
  /*--- Subset the data to keep only observations from the turbine 
        data set. Do not need to keep locations of plants that have
        no wind turbines. */
  if a then do; 
    /*--- Write data set with sites with total capacity exceeding 200 MW. */
    if MW_total > 200 then output US_GT_200MW;
    /*--- Write data set containing all sites in New England. */
    if state in('CT', 'RI', 'VT', 'NH', 'ME', 'MA') then output New_England_All;
  end;
run;

/*--- Generate U.S. map of wind power sites having more than 200 MW capacity.
      Use Open Street Amp background. Specify numeric latitude/longitude
      variables for bubbles and size them based on total wind generation
      capacity at each site (MW_total). */
title 'U.S. Wind Power Sites >200 MW Capacity';
footnote j=l 'U.S. Energy Information Administration (2016)'
         j=r 'Relative MW per site';
proc sgmap plotdata=US_GT_200MW;
  openstreetmap;
  bubble x=long_num y=lat_num size=MW_total;
run;

/*--- Generate map of wind power sites in New England. */
title 'New England Wind Power Sites';
proc sgmap plotdata=New_England_All;
  openstreetmap;
  bubble x=long_num y=lat_num size=MW_total;
run;
quit;