################################################################################ ################################################################################ #Body mass code in R: script to calculate body heights, volume and mass from body length and width measurements of right whales from aerial photographs #Developed by: Fredrik Christiansen (f.christiansen@aias.au.dk, f.christiansen@live.se) #Title: Estimating body mass of free-living whales using aerial photogrammetry and 3D volumetrics #Authors: Fredrik Christiansen, Mariano Sironi, Michael J Moore, Matías Di Martino, Marcos Ricciardi, Hunter A Warick, Duncan J Irschick, Robert Gutierrez & Marcela M Uhart ################################################################################ ################################################################################ ### Import the body morphometric data BM<-data.frame( #Create a data frame with your input data Total.length.m=c(L1,L2,L3,...Ln), #A vector of the body lengths (in meters) of the whales Rep.class=c(R1,R2,R3,...Ln), #A vector of the reproductive classes ("Calf", "Juvenile" or "Lactating" female) of the whales Width.5.proc.m=c(W1,W2,W3...Wn), #A vector of the body widths (in meters) of the whales at 5%BL from the rostrum Width.10.proc.m=c(W1,W2,W3...Wn), #-||- 10%BL -||- Width.15.proc.m=c(W1,W2,W3...Wn), #-||- 15%BL -||- Width.20.proc.m=c(W1,W2,W3...Wn), #-||- 20%BL -||- Width.25.proc.m=c(W1,W2,W3...Wn), #-||- 25%BL -||- Width.30.proc.m=c(W1,W2,W3...Wn), #-||- 30%BL -||- Width.35.proc.m=c(W1,W2,W3...Wn), #-||- 35%BL -||- Width.40.proc.m=c(W1,W2,W3...Wn), #-||- 40%BL -||- Width.45.proc.m=c(W1,W2,W3...Wn), #-||- 45%BL -||- Width.50.proc.m=c(W1,W2,W3...Wn), #-||- 50%BL -||- Width.55.proc.m=c(W1,W2,W3...Wn), #-||- 55%BL -||- Width.60.proc.m=c(W1,W2,W3...Wn), #-||- 60%BL -||- Width.65.proc.m=c(W1,W2,W3...Wn), #-||- 65%BL -||- Width.70.proc.m=c(W1,W2,W3...Wn), #-||- 70%BL -||- Width.75.proc.m=c(W1,W2,W3...Wn), #-||- 75%BL -||- Width.80.proc.m=c(W1,W2,W3...Wn), #-||- 80%BL -||- Width.85.proc.m=c(W1,W2,W3...Wn)) #-||- 85%BL -||- ### Create a table for the height-width ratio for juveniles and lactating females (Table S1) HW.ratios.juvenile.and.lactating<-data.frame( Measurement.site=seq(5,85,5), Juveniles.HW.ratio=c(0.84,0.91,0.96,0.93,0.86,0.86,0.88,0.93,0.99,1.05,1.11,1.15,1.25,1.37,1.56,1.85,2.12), Lactating.HW.ratio=c(0.86,0.90,0.95,0.95,0.83,0.83,0.86,0.92,0.98,1.01,1.04,1.08,1.15,1.23,1.38,1.66,2.10)) ### Create a table for the height-width ratio model parameters for calves (Table S2) HW.ratio.formula.calves<-data.frame( Measurement.site=seq(5,85,5), Intercept=c(1.037,1.113,1.122,1.039,1.198,1.272,1.386,1.575,1.791,2.059,2.350,2.483,2.797,3.176,3.056,2.926,2.206), Slope=c(-0.008,-0.008,-0.009,-0.010,-0.019,-0.022,-0.025,-0.033,-0.042,-0.055,-0.069,-0.073,-0.086,-0.100,-0.085,-0.066,-0.023)) ### Calculate the body heights of whales between 5 and 85%BL from the rostrum for(i in 1:length(BM[,1])){ #Runs a loop for every individual in the data frame temp.widths<-BM[i,c(which(colnames(BM)=="Width.5.proc.m"):which(colnames(BM)=="Width.85.proc.m"))] #Extracts the body width measurements for individual i ifelse(BM$Rep.class[i]=="Calf",{ #If individual i is a calf... Width.55.proc.rel<-BM$Width.55.proc.m[i]/BM$Total.length.m[i]*100 #...extract the relative width (in %) of individual i at 55% body length from the rostrum and... temp.heights<-temp.widths*(HW.ratio.formula.calves$Intercept+HW.ratio.formula.calves$Slope*Width.55.proc.rel) #...calculate the height of calf i from the HW ratio, based on the relative body width of the animal at 55%BL from the rostrum, otherwise... },ifelse(BM$Rep.class[i]=="Juvenile", #If individual i is a juvenile... temp.heights<-temp.widths*HW.ratios.juvenile.and.lactating$Juveniles.HW.ratio, #...calculate the height of juvenile i based on the estimated HW ratios for juveniles, otherwise... temp.heights<-temp.widths*HW.ratios.juvenile.and.lactating$Lactating.HW.ratio)) #If individual i is a lactating female, calculate the height of female i based on the estimated HW ratios for lactating females colnames(temp.heights)<-c("Height.5.proc.m","Height.10.proc.m","Height.15.proc.m","Height.20.proc.m","Height.25.proc.m","Height.30.proc.m","Height.35.proc.m","Height.40.proc.m","Height.45.proc.m","Height.50.proc.m","Height.55.proc.m","Height.60.proc.m","Height.65.proc.m","Height.70.proc.m","Height.75.proc.m","Height.80.proc.m","Height.85.proc.m") #Re-name the column headers for the height measurements ifelse(i==1,temp.output<-temp.heights,temp.output<-rbind(temp.output,temp.heights)) #Extracts the data into a temporary output data frame called "temp.output" } BM<-cbind(BM,temp.output) #Add the predicted height measurements to your main data frame ### Set width and height at 0 and 100%BL from rostrum to be 0 BM$Width.0.proc.m<-0 BM$Height.0.proc.m<-0 BM$Width.100.proc.m<-0 BM$Height.100.proc.m<-0 ### Calculate the body width and height at 90 and 95%BL from rostrum based on linear interpolation between 85 and 100%BL from rostrum BM$Width.90.proc.m<-BM$Width.85.proc.m-(1*(BM$Width.85.proc.m/3)) BM$Width.95.proc.m<-BM$Width.85.proc.m-(2*(BM$Width.85.proc.m/3)) BM$Height.90.proc.m<-BM$Height.85.proc.m-(1*(BM$Height.85.proc.m/3)) BM$Height.95.proc.m<-BM$Height.85.proc.m-(2*(BM$Height.85.proc.m/3)) ### Re-order data frame (necessary for the next step) BM<-BM[,c(which(colnames(BM)=="Width.0.proc.m"),which(colnames(BM)=="Width.5.proc.m"):which(colnames(BM)=="Width.85.proc.m"),which(colnames(BM)=="Width.90.proc.m"),which(colnames(BM)=="Width.95.proc.m"),which(colnames(BM)=="Width.100.proc.m"),which(colnames(BM)=="Height.0.proc.m"),which(colnames(BM)=="Height.5.proc.m"):which(colnames(BM)=="Height.85.proc.m"),which(colnames(BM)=="Height.90.proc.m"),which(colnames(BM)=="Height.95.proc.m"),which(colnames(BM)=="Height.100.proc.m"),which(colnames(BM)=="Total.length.m"),which(colnames(BM)=="Rep.class"))] Width.col.start<-which(colnames(BM)=="Width.0.proc.m") #Extracts the column number of the starting width measurement Height.col.start<-which(colnames(BM)=="Height.0.proc.m") #Extracts the column number of the starting height measurement ### Calculate the body volume of the whales BM$Volume.m3<-NA #Creates an empty storage vector for body volume for(i in 1:length(BM[,1])){ #Runs a loop for every individual in the data frame for(j in 1:(length(seq(0,100,5))-1)){ #Runs a loop for every body segment (volume between two measurement sites), including 0 and 100% f.ellipse<-function(x){ #Formula to calculate the volume of an ellipse (BM[i,Width.col.start+(j-1)]+((BM[i,(Width.col.start+j)]-BM[i,Width.col.start+(j-1)])*x))/2*(BM[i,Height.col.start+(j-1)]+((BM[i,(Height.col.start+j)]-BM[i,Height.col.start+(j-1)])*x))/2*pi } Volume.temp<-integrate(f.ellipse,lower=0,upper=1)$value*BM$Total.length.m[i]*0.05 #Multiplies the area with the volume of the segment ifelse(j==1,Store1<-Volume.temp,Store1<-c(Store1,Volume.temp)) #Stores the output from each body segment volume estimate } BM$Volume.m3[i]<-sum(Store1) #Calculates the body volume from all the body segements for individual i } ### Calculate the body mass of the whales BM$Weight.kg<-BM$Volume.m3*754.6282 #Converts body volume to mass ################################################################################ ################################################################################ #End of script! ################################################################################ ################################################################################