* generate sampling weights for submersed aquatic vegetation;
* created 9/23/09 by BRG;
* last revised: 9/23/09;

* notes:
		- users will want to replace the library name "webcode" used in this code with their own library name
		- for the vegetation component, sampling units correspond to sites rather than to rake surveys. consequently, when datasets contain one record 
				per rake survey (as here), an adjustment must be made to ensure sampling weights are derived from numbers of sites;

%let species = VAAM3; * enter species of interest here;
%let FS = 3;					* enter field station corresponding to trend pool of interest here;
%let libname = webcode;	* enter source and destination location;

* enter population sizes from http://www.umesc.usgs.gov/ltrmp/stats/population_sizes.xls;
data capN_h_VEG;
input fldnum _total_ mstratum $5.;
do year=1998 to 2008 by 1; output; end;	* generate a set of _total_ values for every year;
datalines;
1 7190 BWC-L
1 2391 BWC-U
1 933 BWI
1 986 MCB-L
1 245 MCB-U
1 1351 SC-L
1 897 SC-U
1 2026 TDL-L
1 5478 TDL-U
2 7686 BWC
2 1051 BWI
2 13719 IMP
2 1039 MCB
2 3868 SC
3 11560 BWC
3 1836 BWI
3 14471 IMP
3 2500 MCB
3 2311 SC
4 985 BWC
4 921 BWI
4 488 IMP
4 3914 MCB
4 2968 SC
6 7237 BWC
6 16102 BWI
6 2763 MCB
6 543 SC
run;

* sort population file. restrict to field station of interest;
proc sort data=capN_h_VEG; by fldnum year mstratum; where fldnum=&FS; run;
* sort data file. eliminate non-trend pools;
proc sort data=&libname..&species.FS&FS out=&species.FS&FS; by fldnum year mstratum; where pool notin ("05","12"); run;	
* create sorted data file that contains only one record per barcode/sampling site;
/*proc sort data=&species&FS out=VEGunique nodupkey; by fldnum year mstratum barcode; run;	*/

* calculate sampling weights;
proc means data=&species.FS&FS noprint; var sitedetect; output out=n_h n=n_h; by fldnum year mstratum; run;
data both; merge capN_h_VEG n_h; by fldnum year mstratum;
data weights; set both; if n_h then sweight = _total_ / n_h; else sweight = .; run;

* weights with missing values correspond to sites that were not sampled or, conceivably, to sites that were deemed wholly terrestrial;
proc print data=weights; where sweight=.; run;

* merge with datafile;
data &libname..&species.FS&FS.wt; merge &species.FS&FS weights; by fldnum year mstratum; run;

* confirm weights sum to _totals_ (check log to confirm no records met the mismatch criteria);
/*proc sort data=webcode.SAVP13wt out=SAVP13wtunique nodupkey; by fldnum year mstratum barcode; run;	*/
proc means data=&libname..&species.FS&FS.wt noprint; id _total_; var sweight; output out=wtcheck sum=sumweight; by fldnum year mstratum; run;
proc print data=wtcheck noobs; where round(sumweight,.1) ne _total_ and sumweight; run;