# delimit;
clear;
set memory 8g;
set more off;
version 10;

* #############################################################################;
* set parameters in following section;
* #############################################################################;
* part 1: directories;
* data directory (where downloaded data csv-files are stored);
local dataDir "/media/usbHardDrive/home/ws2162/fineScaleWeather/dataTEXT";
* directory where tempory files are saved;
local tempDir "/tmp";
* directory where end results are saved;
local saveDir "/media/usbHardDrive/home/ws2162/fineScaleWeather/dataSTATA";

* part 2: degree days bound list;
local boundList 0 5 8 10 12 15 20 25 29 30 31 32 33 34;

* part 3: range of years in analysis (can be between 1950 and 2005);
local yearMin = 1950;
local yearMax = 2005;

* part 4: define growing season;
* first month and day of month used in a year;
local monthBegin = 1;  local dayMonthBegin =  1;
* last month and day of month used in a year;
local monthEnd   = 12; local dayMonthEnd   = 31;

* part 5: states that are included in analysis;
local stateList;
* NOTE: delete star at beginning of line to include a state;
* local stateList `stateList'  1; * Arkansas;
* local stateList `stateList'  4; * Arizona;
* local stateList `stateList'  5; * Arkansas;
* local stateList `stateList'  6; * California;
* local stateList `stateList'  8; * Colorado;
* local stateList `stateList'  9; * Connecticut;
* local stateList `stateList' 10; * Delaware;
* local stateList `stateList' 11; * District of Columbia;
* local stateList `stateList' 12; * Folorida;
* local stateList `stateList' 13; * Georgia;
* local stateList `stateList' 16; * Idaho;
* local stateList `stateList' 17; * Illinois;
* local stateList `stateList' 18; * Indiana;
* local stateList `stateList' 19; * Iowa;
* local stateList `stateList' 20; * Kansas;
* local stateList `stateList' 21; * Kentucky;
* local stateList `stateList' 22; * Louisiana;
* local stateList `stateList' 23; * Maine;
* local stateList `stateList' 24; * Maryland;
* local stateList `stateList' 25; * Massachusetts;
* local stateList `stateList' 26; * Michigan;
* local stateList `stateList' 27; * Minnesota;
* local stateList `stateList' 28; * Mississippi;
* local stateList `stateList' 29; * Missouri;
* local stateList `stateList' 30; * Montana;
* local stateList `stateList' 31; * Nebraska;
* local stateList `stateList' 32; * Nevada;
* local stateList `stateList' 33; * New Hampshire;
* local stateList `stateList' 34; * New Jersey;
* local stateList `stateList' 35; * New Mexico;
* local stateList `stateList' 36; * New York;
* local stateList `stateList' 37; * North Carolina;
* local stateList `stateList' 38; * North Dakota;
* local stateList `stateList' 39; * Ohio;
* local stateList `stateList' 40; * Oklahoma;
* local stateList `stateList' 41; * Oregon;
* local stateList `stateList' 42; * Pennsylvania;
* local stateList `stateList' 44; * Rhose Island;
* local stateList `stateList' 45; * South Carolina;
* local stateList `stateList' 46; * South Dakota;
* local stateList `stateList' 47; * Tennessee;
* local stateList `stateList' 48; * Texas;
* local stateList `stateList' 49; * Utah;
* local stateList `stateList' 50; * Vermont;
* local stateList `stateList' 51; * Virginia;
* local stateList `stateList' 53; * Washington;
* local stateList `stateList' 54; * West Virginia;
* local stateList `stateList' 55; * Wisconsin;
* local stateList `stateList' 56; * Wyoming;

* part 6: define the weighting variable and aggregation variable;
local weightVar cropArea;
local aggregationVar fips;

* part 7: define whether you want yearly or monthly output;
* unstar next line for monthly time scale;
local timeScale year month;
* unstar next line for yearly time scale;
* local timeScale year;
* #############################################################################;


* #############################################################################;
* program code from here on, no need to change anything;
* #############################################################################;
* check data - year range;
if (`yearMin' < 1950) {; local yearMin = 1950; };
if (`yearMax' > 2005) {; local yearMax = 2005; };
if (`yearMin' > `yearMax') {;
   local yearMin2 = `yearMin';
   local yearMin = `yearMax';
   local yearMax = `yearMin2';
}; 

******************************************************************************;
* load grid information;
insheet using `dataDir'/gridInfo.csv, clear case;
gen longitude = -125 + mod(gridNumber-1,1405)/24;
gen latitude  = 49.9375+1/48 - ceil(gridNumber/1405)/24;
compress;
sort gridNumber;
save `tempDir'/gridInfo, replace;

******************************************************************************;
* delete possible old files;
! rm `saveDir'/dday_`aggregationVar'_`weightVar'Weighted.dta;

******************************************************************************;
* loop over states;
foreach s of local stateList {;
    * get gridInfo for the state;
    use `tempDir'/gridInfo, clear;
    keep if (floor(fips/1000) == `s');
    keep gridNumber `weightVar' `aggregationVar';
    sort gridNumber;
    save `tempDir'/gridInfo_state, replace;

    * loop over years;
    forvalues t = `yearMin'/`yearMax' {;
        ******************************************************************************;
        * load data;
        insheet using `dataDir'/state`s'_year`t'.csv, clear case;
        gen tAvg = (tMin+tMax)/2;

        * keep only days for selected part of year;
        qui drop if (month < `monthBegin');
        qui drop if (month > `monthEnd');
        qui drop if (month == `monthBegin') & (day < `dayMonthBegin');
        qui drop if (month == `monthEnd')   & (day > `dayMonthEnd');

        * merge in grid weights;
        sort gridNumber;
        merge gridNumber using `tempDir'/gridInfo_state;

        * make sure there is no (_merge == 1) as gridInfo should have data for each gridNumber;
        qui summ _merge;
        if (r(min) == 1) {;
             display("Merge of 1 should not happen - check!")
             break;
        };

        * keep only were weighst are nonzero;
        keep if (_merge == 3) & (`weightVar' > 0);
        drop _merge;

        if (_N > 0) {;
            ******************************************************************************;
            * generate degree days;
            foreach b of local boundList {;
                * create label fro negative bounds by adding Minus;
                if (`b' < 0) {; local b2 = abs(`b'); local bLabel Minus`b2'; }; else {; local bLabel `b'; };

                * default case 1: tMax <= bound;
                gen dday`bLabel'C = 0;

                * case 2: bound <= tMin;
                replace dday`bLabel'C = tAvg - `b' if (`b' <= tMin);

                * case 3: tMin < bound < tMax;
                gen tempSave = acos( (2*`b'-tMax-tMin)/(tMax-tMin) );
                replace dday`bLabel'C = ( (tAvg-`b')*tempSave + (tMax-tMin)*sin(tempSave)/2 )/_pi 
                        if ( (tMin < `b') & (`b' < tMax) );
                drop tempSave;
            };

            ******************************************************************************;
            * collapse by aggregation level;
            collapse tMin tMax tAvg prec dday* [aweight=`weightVar'], by(`aggregationVar' year month day);

            * collapse by year;
            collapse (sum) prec dday* (mean) tMin tMax tAvg, by(`aggregationVar' `timeScale');

            ******************************************************************************;
            * save data;
            capture append using `saveDir'/dday_`aggregationVar'_`weightVar'Weighted;
            sort `aggregationVar' `timeScale';
            save `saveDir'/dday_`aggregationVar'_`weightVar'Weighted, replace;
        };
    };
};

* clear up files;
! rm `tempDir'/gridInfo.dta;
! rm `tempDir'/gridInfo_state.dta;