Stocking Guide in R

[Phil Radtke]

Hi Forest-R:

Anyone interested in taking a look at a script I'm woprking on that
(should) generate stocking guides like the ones developed by Gingrich
(1967) and others (Krajicek et al. 1961, Roach and Gingrich 1968) for
central states (U.S.) hardwoods? My main problem is that using par
(srt=angle) to place labels within the graph only sets angles with
respect to the plot margins, e.g., par(srt=45) rotates a text label
45°. Thus, any rescaling of axes or any change in the parameters of
the stocking guide (especially Dq_min & Dq_max) requires that I tinker
with the rotation angles. The script 'as-is' approximately reproduces
the Gingrich guide for quadratic mean tree diameters between 7 and 15
(Dqmin=7; Dq_max=15) <img src="=http://charcoal.fw.vt.edu/wordpress/wp-
content/uploads/2009/11/StockingGuide.jpeg">, but when I change
Dq_min=3 and Dq_max=7 to match the Gingrich guide for smaller trees,
the limitation of my script becomes evident <img src="http://
charcoal.fw.vt.edu/wordpress/wp-content/uploads/2009/11/
StockingGuide21.jpeg">. I can fix it by tweaking sp_rot, ab_rot, and
Dq_rot, but I'd like to find a way to have the script automatically
set the rotation angles for these labels.

-Phil

# see http://charcoal.fw.vt.edu/wordpress for references
# R-code to generate stocking guide(s) from open-grown crown-width to
DBH relationship
# coefficients of SLR for open grown crown width versus dbh
# from Krajicek (1961) - oak hickory in Central States (For Sci
7:35-42)
a = 3.12;
b = 1.829;

# from Jungkee Choi - Q. varibilis & Q. mongolica
#a = 1.968873754/;
#b = 1.539415729;

# Dq (quadratic mean dbh) limits on the guide size
Dq_max = 15;
Dq_min = 7;

# Rotation angle for stocking percent labels
sp_rot = 42;
# Rotation angle for average tree diameter labels
Dq_rot = -17;
# Rotation angle for A and B labels
ab_rot = -5;

# define functions
tpa <- function(ccf,dq) (ccf*43560)/(100*((a+b*dq)/2)^2*pi)
ba <- function(ccf,tpa)
{
((sqrt((ccf*43560)/(100*tpa*pi/4))-a)/b)^2*tpa*pi/(4*144);
}

# generate upper limit for scaling the stocking guide axes
CCF_max = 200;
CCF_min = 70;
CCF_lower = 87.5;
CCF_upper = 192.5;
CCF_limits <- c(CCF_lower,CCF_upper);
TPA_max = tpa(CCF_max,Dq_min);
TPA_min = tpa(CCF_min,Dq_max);
BA_max <- ba(CCF_max,tpa(CCF_max,Dq_max));
BA_min <- ba(CCF_min,tpa(CCF_min,Dq_min));

# generate data for the stocking guide B line
CCF_b_line = 100;
TPA_max_b_line = tpa(CCF_b_line,Dq_min);
TPA_min_b_line = tpa(CCF_b_line,Dq_max);
TPA_b_line <- seq(TPA_min_b_line,TPA_max_b_line,length=100);
BA_b_line <- ba(CCF_b_line,TPA_b_line);

# generate data for the stocking guide A line
CCF_a_line = 175;
TPA_max_a_line = tpa(CCF_a_line,Dq_min);
TPA_min_a_line = tpa(CCF_a_line,Dq_max);
TPA_a_line <- seq(TPA_min_a_line,TPA_max_a_line,length=100);
BA_a_line <- ba(CCF_a_line,TPA_a_line);

# plot A line and B line in a scaled X-Y window
plot(TPA_a_line,BA_a_line,type='l',
xlab="Trees per acre",
ylab="Basal area per acre (square feet)",
xlim=c(TPA_min,TPA_max),
ylim=c(BA_min,BA_max),lwd=3);
par(xaxp=c(par("xaxp")[1],par("xaxp")[2],par("xaxp")[3]*2));
par(yaxp=c(par("yaxp")[1],par("yaxp")[2],par("yaxp")[3]*2));
axis(1);
axis(2);
par(xaxp=c(par("xaxp")[1],par("xaxp")[2],par("xaxp")[3]*2));
par(yaxp=c(par("yaxp")[1],par("yaxp")[2],par("yaxp")[3]*2));

# shade between A and B lines light gray
polygon(c(TPA_a_line,sort(TPA_b_line,decreasing=TRUE)),
c(BA_a_line,sort(BA_b_line,decreasing=FALSE)),col=gray(.
9),bg="transparent");
grid();

# shade above A line medium gray
TPA_max = tpa(CCF_upper,Dq_min);
TPA_min = tpa(CCF_upper,Dq_max);
TPA <- seq(TPA_min,TPA_max,length=100);
BA <- ba(CCF_upper,TPA);
polygon(c(TPA_a_line,sort(TPA,decreasing=TRUE)),
c(BA_a_line,sort(BA,decreasing=FALSE)),col=gray(.
6),bg="transparent");
grid();
lines(TPA_b_line,BA_b_line,lwd=3);
lines(TPA_a_line,BA_a_line,lwd=3);

# label A and B lines
par(srt=ab_rot);
text(min(TPA_b_line),max(BA_b_line),labels="B",cex=1.1,font=2,adj=c
(1.2,.4));
text(min(TPA_a_line),max(BA_a_line),labels="A",cex=1.1,font=2,adj=c
(1.2,.4));

# draws stocking percent lines at intervals of 17.5% of CCF
# à la Gingrich (1967, For Sci 13:38-53)
for(i in seq(CCF_lower,CCF_upper,length=7)) {
TPA_max = tpa(i,Dq_min);
TPA_min = tpa(i,Dq_max);
TPA <- seq(TPA_min,TPA_max,length=100);
BA <- ba(i,TPA);
lines(TPA,BA,lwd=1);
par(srt=sp_rot); # rotating this at a fixed angle seems sub-optimal
text(TPA_max,ba(i,TPA_max),labels=(i/1.75),adj=c(.5,1.2),cex=0.9)
if(i == CCF_lower + 17.5)
{
par(srt=sp_rot);
text(TPA_max,ba(i,TPA_max),labels="STOCKING PERCENT",adj=c
(0,2.5))
}
}
# draws the lines representing 1" intervals of quadrratic mean
diameter
for(i in seq(Dq_min,Dq_max,length=(Dq_max-Dq_min+1))) {
TPA <- tpa(CCF_limits,i);
BA <- ba(CCF_limits,TPA); lines(TPA,BA,lwd=1); par
(srt=Dq_rot); # rotating this at a fixed angle seems sub-optimal
text(TPA[2],BA[2],labels=(i),adj=c(.6,-.3),cex=0.9) if(i > ((Dq_max
+Dq_min)/2) && i <= (((Dq_max+Dq_min)/2)+1))
{
par(srt=Dq_rot);
text(TPA[2],BA[2]*1.05,labels="AVERAGE TREE DIAMETER",adj=0);
text(sum(TPA)/1.3,sum(BA)/1.57,labels="OVERSTOCKED",adj=0);
text(sum(TPA)/1.8,sum(BA)/1.95,labels="FULLY STOCKED",adj=0);
text(TPA[1],BA[1]*1.1,labels="UNDERSTOCKED",adj=0);
}
}