* Redirect SAS output from output window to PDF file;
ods listing close;
ods noptitle;
ods pdf file='D:\SURVEY\Pws2013\Western PWS\Western PWS 2013.pdf' notoc;

%let st=sightings_identity_rep1_merge;									/* set sighting table */
%let mdb=D:\SURVEY\Pws2013\Western PWS\Western_PWS_2013.mdb;			/* set geodatabase */
%let tt=transects_identity_rep1_merge;									/* set transect table */
%let h=1003.35; 														/* set high density stratum size */
%let l=1354.51; 														/* set low density stratum size  */

options nodate pagesize=66 linesize=80 nonumber orientation=portrait;
*libname test 'C:\Esslinger\SURVEY\Pws2012\Western PWS\Western_PWS_2013.mdb'; 	/* geodatabase from ArcMap */

title1 'US Geological Survey';
title3 'Alaska Science Center';
title5 '2013 Western PWS Sea Otter Survey';

***************************************************************************************;
*   Sighting data;
***************************************************************************************;
data one;
proc import dbms=accesscs
			out=test replace
			table="&st";												/* import sightings within REPLICATE survey area */
			database="&mdb";
run;

* Subset strip count data into regular and complete counts;
data sighting complete allspecies; 										/* define 3 output data sets */
  set test; 															/* read from Access table in geodatabase */
    keep transect stratum obs_datetime species group  					/* keep these fields for analysis */
         strip_adults strip_pups strip_count 
         circle_adults circle_pups circle_count
         FID_Sightings;

	  array sa{10} SA01-SA10; 											/* define array for strip adult fields */
	  array sp{10} SP01-SP10; 											/* define array for strip pup fields */
	  array ca{10} CA01-CA10; 											/* define array for circle adult fields */
	  array cp{10} CP01-CP10; 											/* define array for circle pup fields */
	  array cc{10} CC01-CC10; 											/* define array for complete count fields */

	  do i=1 to 10;														/* initialize DO loop from 1-10 */

	    group=i; 														/* set group ID value to i */

		if sa{i}=-999 then sa{i}=0;										/* if strip adult field is null, set to zero */
		if sp{i}=-999 then sp{i}=0;										/* if strip pup field is null, set to zero */
		if ca{i}=-999 then ca{i}=.;										/* if circle adult field is null, set to missing */
		if cp{i}=-999 then cp{i}=.;										/* if circle pup field is null, set to missing */

	    strip_adults=sa{i};												/* set strip_adult value for group i */
		strip_pups=sp{i}; 												/* set strip_pup value for group i */
		strip_count=strip_adults+strip_pups; 							/* calculate strip_count value for group i */

		circle_adults=ca{i}; 											/* set circle_adult value for group i */
		circle_pups=cp{i}; 												/* set circle_pup value for group i */
		circle_count=circle_adults+circle_pups;							/* calculate circle_count value for group i */

		if cc{i}=0 and species IN("SEOT","ISU") then do; 				/* if not a complete count of sea otters, then do */				
		   if strip_count > 0 or circle_count > 0  						/* if strip_count or circle_count greater than zero */
              then output sighting;										/* then output record to SIGHTING dataset */
		end;

		if cc{i}=1 and species IN("SEOT","ISU") then do; 				/* if it is a complete count of sea otters then do */						
           if strip_count > 0 or circle_count > 0  						/* if strip_count or circle_count >greater than zero */
              then output complete;										/* then output record to COMPLETE dataset */
		end;

		if strip_count > 0 then output allspecies; 						/* if strip_count greater than zero, then output to ALLSPECIES dataset */

	  end;
run;

proc sort data=sighting; 												/* sort SIGHTING dataset by observation datetime */
  by obs_datetime;
run;

proc sort data=complete; 												/* sort SIGHTING dataset by observation datetime */
  by obs_datetime;
run;

data strip_sighting; 													/* create STRIP_SIGHTING dataset for strip_count values greater than zero */
  set sighting;
    if strip_count=0 then delete;
run;

* Total otters by transect;
proc summary sum nway data=strip_sighting; 		
  class transect;
    var strip_count;
	    id stratum;
          output out=sightsum sum=;										/* output to SIGHTSUM dataset */
run; 
proc summary sum nway data=complete;
  class transect;
    var strip_count;
	    id stratum;
          output out=compsum sum=; 										/* output to COMPSUM dataset */
run;

***************************************************************************************;
*   Transect data;
***************************************************************************************;
data two;
proc import dbms=accesscs
			out=test2 replace
			table="&tt";												/* import transects within REPLICATE survey area */
			database="&mdb";	
run;

data transect; 
  set test2; 															/* read from Access table in geodatabase */
    keep transect stratum survey_date length area; 						/* keep these fields for analysis */
	  length=shape_length; 												/* create new LENGTH field (simpler to remember) */
	  if survey_date = . then delete; 									/* drop unsurveyed transects   */
      area=(length*.0004);     											/* calculate area in square km */
run;

proc sort data=transect;  												/* sort by transect */
  by transect;
run;

***************************************************************************************;
*   Combine Sighting and Transect data;
***************************************************************************************;
data sightsum;
  merge sightsum (in=insight) transect;  								/* merge SIGHTING and TRANSECT datasets */
    by transect;  														/* by transect field */
	  drop _type_; 														/* drop _TYPE_ field, it's annoying */
	    if strip_count=. then strip_count=0; 							/* if strip_count is missing, set to zero for transect */
		if _freq_ = . then _freq_ = 0; 									/* if _freq_ is missing, set to zero for transect */
		label stratum='Survey Stratum'									/* label variables for nicer output */
              length='Transect Length (m)'
              area='Transect Area (km^2)'
              _freq_='Otter Groups (n)'
              strip_count='Total Otters';
run;

data compsum; 															/* do the same as above to merge COMPSUM and TRANSECT datasets */
  merge compsum (in=incomp) transect;
    by transect;
	  drop _type_;
	    if strip_count=. then strip_count=0;
		if _freq_ = . then _freq_ = 0;
	    label stratum='Survey Stratum'
              length='Transect Length (m)'
              area='Transect Area (km^2)'
              _freq_='Otter Groups (n)'
              strip_count='Total Otters';
run;

proc sort data=sightsum; by stratum; run;

title9 'Transect and Small Group Sighting Information';
proc print noobs label data=sightsum; 									/* print strip sighting information for all transects */
  by stratum; 															/* separate by stratum */
     var transect _freq_ strip_count length area; 						/* include these fields */
	    sum strip_count length area; 									/* summarize these fields */
run;
title9 'Transect and Complete Count Information';
proc print noobs label data=compsum; 									/* print complete count information */
  by stratum;
    var transect _freq_ strip_count length area;
	  sum strip_count length area;
	    where strip_count > 0; 											/* only print transects that had complete counts (saves paper) */
run;

***************************************************************************************;
*   Calculate ratio estimates of density;
***************************************************************************************;

proc summary nway sum data=sightsum; 									/* summarize strip data */
  class stratum; 														/* by stratum */
    var strip_count area; 												/* summarize strip_count and area */
      output out=density1 sum=totcount totarea;							/* output to DENSITY1 dataset */
run;

proc summary nway sum data=compsum;										/* do the same as above for complete count data */
  class stratum;
    var strip_count area;
      output out=density2 sum=totcount totarea;
run;

data density1; 													
  set density1;															/* process strip count data */
    if totcount > 0 then  												/* if total count is greater than zero */
       density=totcount/totarea; 										/* then calculate density in otters/km^2 */
	else density = 0; 													/* otherwise, density equals zero */
	  drop _type_; 														/* drop this annoying _type_ variable again.  I mean, really? */
	  rename _freq_ = num; 												/* rename _freq_ to num just because */
run;

data density2; 															/* do the same for complete count data */
  set density2;
    if totcount > 0 then 
       density=totcount/totarea;
	else density = 0;
	  drop _type_;
      rename _freq_ = num;
run;

data sightsum2;															
  merge sightsum density1; 												/* merge strip count data with density value */
    by stratum;															/* by stratum */

      sumsqr+((strip_count-area*density)**2); 							/* increment sum-of-squared errors */

	    if stratum='H' then strarea=&h;									/* set high density stratum area with macro variable */
		if stratum='L' then strarea=&l;									/* set low density stratum area with macro variable */

		if last.stratum = 1 then do; 									/* if it's the last record in the stratum then do */
		  popsize=density*strarea; 										/* population size equals density x stratum area */
		  fract=totarea/strarea; 										/* fraction surveyed equals sampled area / total area */
          pop_var=(strarea**2)*(1-fract)*sumsqr/((totarea**2/num)*(num-1));		/* population variance equals this thing */
          stderr=pop_var**.5; 											/* standard error equals square root of variance */
          pstderr=stderr/popsize; 										/* coefficient of variation equals standard error / population size */
		  output; 														/* output this last record to the dataset */
		  sumsqr=0;														/* reset sum-of-squared errors to zero */
		end; 															/* WHEW! */

		label totarea='Total Area (km^2)'								/* label these variables or you won't know what the heck they mean */
		      totcount='Total Count'
			  num='Transects (n)'
			  density='Density (otters/km^2)'
			  sumsqr='Sum of Squared Error'
			  strarea='Stratum Area'
			  stratum='Survey Stratum'
			  popsize='Uncorrected Population Size'
			  fract='Percent Sampled'
			  pop_var='Population Variance'
			  stderr='Population Standard Error'
			  pstderr='Coefficient of Variation';
run;

title9 'Uncorrected Population Estimates';

proc print data=sightsum2 noobs label; 									/* print uncorrected estimates by stratum */
  var stratum num totcount totarea density strarea popsize pop_var stderr pstderr;
    where stratum <> '';
run;

proc sort data=compsum; by stratum; run;

data compsum2; 															/* do the same as above for complete count data */
  merge compsum density2; 												/* you don't really need comments here, right? */
    by stratum;															/* i mean, just look up a little */

      sumsqr+((strip_count-area*density)**2);

	    if stratum='H' then strarea=&h;
		if stratum='L' then strarea=&l;

		if last.stratum = 1 then do;
		  popsize=density*strarea;
		  fract=totarea/strarea;
		  pop_var=(strarea**2)*(1-fract)*sumsqr/((totarea**2/num)*(num-1));
		  stderr=pop_var**.5;
          if popsize > 0 then pstderr=stderr/popsize;
		  output;
		  sumsqr=0;
		end; 															/* hope you're not mad about the comments... */

		label totarea='Total Area (km^2)'
		      totcount='Total Count'
			  num='N'
			  density='Density (otters/km^2)'
			  sumsqr='Sum of Squared Error'
			  strarea='Stratum Area'
			  stratum='Survey Stratum'
			  popsize='Population Size'
			  fract='Percent Sampled'
			  pop_var='Population Variance'
			  stderr='Population Standard Error'
			  pstderr='Coefficient of Variation';
run;

title9 'Complete Count Population Estimates';

proc print data=compsum2 noobs label; 									/* print complete count estimates by stratum */
  var stratum totcount totarea density strarea popsize pop_var stderr pstderr;
     where stratum <> '';
run;


***************************************************************************************;
*   ISU data;
***************************************************************************************;
data isu1; 																/* create new ISU1 dataset from SIGHTING dataset */
  set sighting;
    retain i; 															/* retain this value for ISU number */
    if _n_=1 then i=0; 													/* initialize ISU number to zero at the first record */

    by obs_datetime; 													/* remember, obs_datetime is a unique identifier for all groups at each sighting location */
	  if first.obs_datetime=1 and circle_count > 0 then do;				/* if the first group at that sighting location, then do */
         i=i+1;															/* increment the ISU number by 1 */
	  end; 													

	if circle_count > 0 then isu=i;										/* if circle_count greater than zerio, then set ISU value */

	drop i;																/* drop i, because that stuff adds up, man! */
run;

title9 'All ISU Information';
proc print noobs label data=isu1; 										/* print all the ISU information to see what you're working with */
  var isu transect group strip_count circle_count;
    where isu <> .;														/* only print records form ISUs */
run;

data isu2; 																/* make ISU2 dataset */
  set isu1;
    if group=1 then delete;												/* drop first group from calculations because Mark Udevitz says so */
run;

proc summary sum nway data=isu2; 										/* summarize ISU data */
  class isu; 															/* by ISU number */
    var strip_count circle_count; 										/* sum strip_count and circle_count variables */
	  id transect;														/* include transect id in output data set */
      output out=isusum sum=; 											/* create ISUSUM dataset */
run;

title9 'Summary of Usable ISU Information';

proc print noobs label data=isusum;										/* print summary of usable ISUs */
   var isu transect strip_count circle_count;
	  sum strip_count circle_count;
run;

* Calculate ratio estimate of correction factor;

proc summary sum nway data=isusum; 										/* summarize all ISU data to calculate correction factor */
   var strip_count circle_count;
      output out=corrfact n=num sum=strip_total circle_total; 			/* output to CORRFACT dataset and rename summary variables */
run;

data corrfact; 						
  set corrfact;
    factor=circle_total/strip_total; 									/* calculate correction factor as circle_total / strip_total */

	do i = 1 to num; 													/* create multiple records for next step */
	  output;
	end;
run;

data corrfact2;  														
  merge isusum corrfact; 												/* merge ISUSUM and CORRFACT datasets */
run;

data corrfact2;						
  set corrfact2;

    sumsqr+((circle_count-strip_count*factor)**2); 						/* increment sum-of-squared errors */
	sumsp+strip_count; 													/* increment strip_count */

	if _n_ = num then do; 												/* if the last record in the dataset then do this stuff */
       fact_var=sumsqr/((sumsp**2/num)*(num-1)); 						/* correction factor variance equals this thing */
       stderr=sqrt(fact_var); 											/* standard error equals square root of variance (Statistics 101, Bro) */
       pstderr=stderr/factor; 											/* coefficient of variation */
	   output; 															/* output this record to the dataset */
	end;

	label num='N ISUs'													/* label these variables because this dataset is just as bad as the other ones */
	      strip_total='Total Strip Count'
		  circle_total='Total Circle Count'
		  factor='Correction Factor'
		  sumsqr='Sum of Squared Error'
		  fact_var='Variance of Correction Factor'
		  stderr='Standard Error of Correction Factor'
		  pstderr='Coefficient of Variation';
run;

title9 'ISU Correction Factor';

proc print noobs label data=corrfact2; 									/* print this correction factor mumbo-jumbo and don't ask too many questions */
   var num strip_total circle_total factor fact_var stderr pstderr;
run;

/******************************************************************/
/* Code to obtain adjusted population size estimates.             */
/******************************************************************/

* Duplicate correction factor for each stratum;

data corrfact3;				
  set corrfact2;
    output;																/* this one is for the high density stratum */
    output;																/* and this one is for the low density stratum */
run; 																	/* or was it the other way around? */

* Calculate adjusted population size estimate, variance, standard error, and proportional standard error; 
data correst; 
  merge sightsum2 corrfact3; 											/* merge the uncorrected estimates with correction factor */
    adjpopsz=popsize*factor; 											/* calculate adjusted population size */
    variance=(popsize**2)*fact_var + (factor**2)*pop_var- pop_var*fact_var;		/* calculate variance of adjusted population size as a product (Statistics 201) */
    stderr=variance**.5; 												/* calculate standard error */
    pstderr=stderr/adjpopsz; 											/* coefficient of variation */

	label adjpopsz='Adjusted Population Size'							/* variable labels */
	      variance='Adjusted Population Variance'
		  stderr='Adjusted Population Standard Error'
		  pstderr='Coefficient of Variation';
run;

title9 'Corrected Population Estimates by Stratum';

proc print noobs label data=correst; 									/* print corrected population estimates by stratum */
  var stratum popsize factor adjpopsz variance stderr pstderr;
    where stratum <> '';
run;

* Code to sum corrected and uncorrected estimates;

data correst; 								
  set correst;
    popsize=adjpopsz;													/* copy adjusted population size variable to add with complete count variable */
	pop_var=variance;													/* copy adjusted variance variable to add with complete count variance */
      keep stratum popsize pop_var;
run;

data complete; 															/* minimize complete count data set to contain only the variables needed */
  set compsum2;
	  keep stratum popsize pop_var;										/* that would be these three */
run;

data correst2; 															/* combine the two data sets */
  set correst complete;
    if popsize=0 then delete; 											/* throw out any records that have population size of zero (usually complete counts for low density stratum */
run;

proc summary sum nway data=correst2; 								
  var popsize pop_var; 													/* add population size and variance together */
    output out=finest sum=;
run;

data finest;
  set finest;
    stderr=pop_var**.5;													/* calculate standard error just one more time */
	pstderr=stderr/popsize;												/* calculate coefficient of variation, too! */

	label stderr='Standard Error'										/* label these variables */
	      pstderr='Coefficient of Variation'
          popsize='Population Estimate'
          variance='Population Variance';
run;

title9 'Final Population Estimate';
proc print noobs label data=finest;										/* print your final population estimate */
  var popsize pop_var stderr pstderr;
run;

/******************************************************************/
/* Make Sighting and ISU table in Geodatabase                     */
/******************************************************************/

proc datasets library=test nolist; 										/* delete existing table from the Access database */
  delete SightSum; 														/* otherwise you'd get an error in your log file */
run; 																	/* nobody likes that */

proc summary data=allspecies n sum nway; 								/* summarize all species data - you were wondering when we'd get to that, weren't you? */
  class FID_Sightings; 													/* add stuff up for each sighting location */
    var strip_count circle_count; 										/* summarize strip_count and circle_count variables */
	  id species;
	  output out=allsum	n= strip_groups circle_groups  					/* count the number of strip and circle groups */
             sum=strip_total circle_total; 								/* add up the total number of strip and circle critters */
run; 	

proc summary data=isu1 n nway; 											/* summarize the ISU information at each sighting location */
  class FID_Sightings; 					
    var isu;
	  id isu;
	    where isu <> .;
        output out=isu1sum n=isugroups;
run;

data test.SightSum;														/* create Access table in geodatabase */
  merge allsum isu1sum;													/* merge all species summary and ISU smmary */
    by FID_Sightings; 													/* by observation date and time */
run;

* Now you will need to open the database in Microsoft Access and create an update query to plug these summary values into the attribute table for your SIGHTINGS feature class;
* You can also delete all those array fields, as you won't be using them for anything in ArcMap;
* In ArcMap, select just those sightings that have an ISU number greater than zero, and buffer them by 200 meters to create an ISU feature for mapping;
* If you really want to get fancy, you can move those ISU circles to line up better with the survey swath - that's up to you;

ods pdf close;															/* stop directing output to the PDF file */
ods listing; 															/* and send it to the LIST window instead */