Skip to content

Commit

Permalink
Merge pull request #9 from mle2718/update2024
Browse files Browse the repository at this point in the history
Update2024
  • Loading branch information
mle2718 authored Apr 22, 2024
2 parents ea2bab4 + 104948e commit 3206417
Show file tree
Hide file tree
Showing 22 changed files with 1,711 additions and 73 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,18 @@
/* read in and store the MRIP location stuff. Put it in the file ma_site_allocation.dta */
clear
#delimit;
odbc load, exec("select site_id, stock_region_calc from recdbs.mrip_ma_site_list;") $myNEFSC_USERS_conn;
#delimit cr
save "${data_raw}/ma_site_allocation.dta", replace
odbc load, exec("select site_id, stock_region_calc from recdbs.mrip_ma_site_list;") $mysole_conn;

/* will need to switch the connection string to $myNEFSC_USERS_conn soon */


/* patch in the extra site_id */

qui count;
local newobs=`r(N)'+1;
set obs `newobs';

replace site_id=4434 if site_id==.;
replace stock_region_calc="NORTH" if site_id==4434;

save "${data_raw}/ma_site_allocation.dta", replace;
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
/* after running do $BLAST_MRIP */
/* get the mass location data from oracle */

*do ${extraction_code}/get_ma_allocation.do
do ${extraction_code}/get_ma_allocation.do


/* copying over partial years of data */
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,243 @@
/* This is a file that produces a dataset that contains #of fish encountered per trip.
This is a port of Scott's "cod_domain_length_freqs_by_wave_gom_2013.sas"
This is a template program for estimating length frequencies
using the MRIP public-use datasets.
The program is setup to use information in the trip_yyyyw
dataset to define custom domains. The length frequencies are
estimated within the domains by merging the trip information onto
the size_yyyyw datasets.
Required input datasets:
trip_yyyyw
size_yyyyw
It looks like we also need to port cod_domain_length_freqs_b2_by_wave_gom_2013 as well
There will be one output per variable and year in working directory:
"$my_common`myv'_b2_OpenClose_$working_year.dta"
*/

/* General strategy
COMPUTE totals and std deviations for cod catch
*/
clear
mata: mata clear




tempfile tl1 sl1

foreach file in $triplist{
append using ${data_raw}/`file'
}
capture drop $drop_conditional


replace var_id=strat_id if strmatch(var_id,"")
sort year strat_id psu_id id_code
/* Deal with new variable names in the transition period
*/

capture confirm variable wp_int_chts
if _rc==0{
drop wp_int
rename wp_int_chts wp_int
else{
}
}
capture confirm variable wp_size_chts
if _rc==0{
drop wp_size
rename wp_size_chts wp_size
else{
}
}
save `tl1'
clear


foreach file in $b2list{
append using ${data_raw}/`file'
}
cap drop $drop_conditional
/* Deal with new variable names in the transition period */

capture confirm variable wp_int_chts
if _rc==0{
drop wp_int
rename wp_int_chts wp_int
else{
}
}
capture confirm variable wp_size_chts
if _rc==0{
drop wp_size
rename wp_size_chts wp_size
else{
}

}

capture confirm variable wp_catch_chts
if _rc==0{
drop wp_catch
rename wp_catch_chts wp_catch
else{
}

}

sort year strat_id psu_id id_code
replace common=subinstr(lower(common)," ","",.)
save `sl1'

use `tl1'
merge 1:m year strat_id psu_id id_code using `sl1', keep(1 3)
drop _merge
keep if year==$working_year
/* THIS IS THE END OF THE DATA MERGING CODE */


/* ensure that domain is sub_reg=4 (New England), relevant states (MA, NH, ME), mode_fx =123, 457 */
keep if sub_reg==4
keep if st==23 | st==33 |st==25

/*This is the "full" mrip data */
tempfile tc1
save `tc1'

/*classify as GOM or GB based on the ma_site_allocation.dta file */
rename intsite site_id
sort site_id


merge m:1 site_id using "${data_raw}/ma_site_allocation.dta", keepusing(stock_region_calc)
rename site_id intsite
drop if _merge==2
drop _merge

/*classify into GOM or GBS */
gen str3 area_s="OTH"

replace area_s="GOM" if st==23 | st==33
replace area_s="GOM" if st==25 & strmatch(stock_region_calc,"NORTH")
replace area_s="GBS" if st==25 & strmatch(stock_region_calc,"SOUTH")

/* classify catch into the things I care about (common==$mycommon) and things I don't care about "ZZZZZZZZ" */
gen common_dom="Z"
if strmatch("$my_common","atlanticcod")==1{
replace common_dom="C" if strmatch(sp_code,"8791030402")
}

if strmatch("$my_common","haddock")==1{
replace common_dom="H" if strmatch(sp_code,"8791031301")
}




tostring wave, gen(w2)
tostring year, gen(year2)

destring month, gen(mymo)
drop month
tostring mymo, gen(month)
drop mymo

tostring year, gen(myy)




/* l_in_bin already defined
gen l_in_bin=floor(lngth*0.03937) */

/* this might speed things up if I re-classify all length=0 for the species I don't care about */
replace l_in_bin=0 if strmatch(common_dom, "Z")==1


sort year w2 strat_id psu_id id_code


/* use unweighted for atlantic cod lengths */
if "$my_common"=="atlanticcod"{
svyset psu_id, strata(var_id) singleunit(certainty)
gen open=inlist(month,"9","10")

}
if "$my_common"=="haddock"{
svyset psu_id [pweight= wp_size], strata(var_id) singleunit(certainty)
gen open=!inlist(month,"3")

}

tostring open, gen(myo)

gen my_dom_id_string=common_dom+"_"+area_s+"_"+myo



local myv l_in_bin

svy: tab `myv' my_dom_id_string, count
/*save some stuff
matrix of proportions, row names, column names, estimate of total population size*/
mat eP=e(Prop)
mat eR=e(Row)'
mat eC=e(Col)
local PopN=e(N_pop)

local mycolnames: colnames(eC)
mat colnames eP=`mycolnames'

clear
/*read the eP into a dataset and convert proportion of population into numbers*/
svmat eP, names(col)
foreach var of varlist *{
replace `var'=`var'*`PopN'
}
/*read in the "row" */
svmat eR
order eR
rename eR `myv'



keep `myv' *GOM*
drop *Z_*
gen year=$working_year


qui desc

foreach var of varlist *GOM*{
tokenize `var', parse("_")
rename `var' `3'_`5'
}
reshape long GOM_, i(year l_in_bin) j(oo)
rename GOM_ count
gen str6 open="OPEN" if oo==1
replace open="CLOSED" if oo==0

sort year open l_in_bin
order year open l_in_bin
drop oo


save "$my_outputdir/$my_common`myv'_b2_OpenClose_${working_year}.dta", replace







Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,13 @@ replace l_in_bin=0 if strmatch(common_dom, "Z")==1
sort year w2 strat_id psu_id id_code
svyset psu_id [pweight= wp_size], strata(var_id) singleunit(certainty)


/* use unweighted for atlantic cod lengths */
if "$my_common"=="atlanticcod"{
svyset psu_id, strata(var_id) singleunit(certainty)
}



local myv l_in_bin

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -201,10 +201,12 @@ svyset psu_id [pweight= wp_int], strata(var_id) singleunit(certainty)
preserve
sort my_dom_id year strat_id psu_id id_code
local myvariables dtrip
local i=1
local i=0
/* total with over(<overvar>) requires a numeric variable */

foreach var of local myvariables{
local ++i

svy: total `var', over(my_dom_id)

mat b`i'=e(b)'
Expand All @@ -213,9 +215,7 @@ foreach var of local myvariables{
mat sub`i'=vecdiag(V)'
mat colnames sub`i'=variance

local ++i
}
local --i
sort my_dom_id year
duplicates drop my_dom_id, force
keep my_dom_id year area_s month dom_id
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@


global cod_lwa 0.000005132
global cod_lwb 3.1625
global had_lwa 0.000009298
global had_lwe 3.0205
global lngcat_offset_cod 0.5
global lngcat_offset_haddock 0.5
global mt_to_kilo=1000
global kilo_to_lbs=2.20462262
global cm_to_inch=0.39370787
global mortality_release=0.15





use "$annual/cod_size_class_OpenClose_$this_working_year.dta" if open=="OPEN", clear
expand 2 , gen(mark)
gen month=9
replace month=10 if mark==1

drop mark
tempfile open
save `open'


use "$annual/cod_size_class_OpenClose_$this_working_year.dta" if open=="CLOSED", clear
expand 12
bysort year open lngcat: gen month=_n
drop if inlist(month,9,10)

append using `open'

drop count
bysort year month open: egen tab1=total(ab1_count)
bysort year month open: egen tb2=total(b2_count)

gen prob_ab1=ab1_count/tab1
gen prob_b2=b2_count/tb2


keep year month open lngcat prob*

tempfile lengths
save `lengths', replace


use "${monthlydir}\atlanticcod_catch_$this_working_year.dta", clear
gen ab1=claim+harvest
rename release b2

gen str6 open="OPEN" if inlist(month,"09","10")
replace open="CLOSED" if open==""
destring month, replace
keep year month ab1 b2 open
merge 1:m year month open using `lengths', keep(1 3)
gen ab1_count=ab1*prob_ab1
gen b2_count=b2*prob_b2
keep year month lngcat open ab1_count b2_count
sort year month open lngcat


gen weight_per_fish= $kilo_to_lbs*$cod_lwa*((lngcat+.5)/$cm_to_inch)^$cod_lwb
gen ab1mt=weight_per_fish*ab1_count/($mt_to_kilo*$kilo_to_lbs)
gen b2mt=weight_per_fish*b2_count/($mt_to_kilo*$kilo_to_lbs)



collapse (sum) ab1_count b2_count ab1mt b2mt , by(year month)

gen b2dead_mt=b2mt*$mortality_release

format ab1_count b2_count %10.0fc

format b2mt ab1mt b2dead_mt %6.2fc

save "$annual\atlanticcod_weights_OpenClose_${this_working_year}.dta", replace


Loading

0 comments on commit 3206417

Please sign in to comment.