library(ggplot2) library(reshape2) ###############PROBLEM WITH facet_wrap including variable names freq_hz makes facet_wrap put in empty squares library(directlabels) library(scales) library(plyr) ExportPlot <- function(gplot, filename, width=2, height=1.5) { # Export plot in PDF and EPS. # Notice that A4: width=11.69, height=8.27 ggsave(paste(filename, '.pdf', sep=""), gplot, width = width, height = height) ggsave(paste(filename, '.jpg', sep=""), gplot, width = width, height = height) postscript(file = paste(filename, '.eps', sep=""), width = width, height = height) print(gplot) dev.off() png(file = paste(filename, '_.png', sep=""), width = width * 100, height = height * 100) print(gplot) dev.off() } ########################## data files load(file="./data_1_3_BB_100.Rdata") # individual mmsi transits head(data_1_3_BB_100) datetime mon hr mmsi shipname classA classB length width deadweight year rnm sog fhz dBrl ddBsig ddBback 1 3/7/2011 18:33:19 3 18 316011773 GRIZZLY 60 Military Military 0 0 0 0 1.16 9.5 0 110 109.92 89.80 24 3/8/2011 1:38:37 3 1 247279500 AURORA D Bulk carrier Bulk carrier 200- 158 26 23923 1997 1.49 14.0 0 110 109.85 93.99 54 3/8/2011 4:52:7 3 4 245855000 MAKIRI GREEN Cargo Cargo 143 22 17539 1999 1.04 14.3 0 114 113.91 96.50 91 3/8/2011 9:4:8 3 9 316005594 OCEAN CLIPPER Tug Tug 27 9 0 0 0.91 3.6 0 100 91.30 98.00 101 3/8/2011 15:43:27 3 15 372359000 IYO WIND Bulk carrier Bulk carrier 200- 190 32 53569 2008 1.33 14.2 0 114 113.94 93.39 139 3/8/2011 20:46:7 3 20 636091891 VANCOUVER EXPRESS Container ship Container ship 320+ 334 43 103773 2009 1.26 22.1 0 117 116.96 94.49 ddB_SL_20 ddB_SL_20_abs ddB_SL_Emp ddB_SL_EmpAbs ddB_SL_JASCO 1 176.56 176.64 171.93 172.00 NaN 24 178.67 184.28 173.88 179.50 NaN 54 179.60 180.36 175.03 175.80 NaN 91 155.84 155.88 151.35 151.39 NaN 101 181.77 196.54 177.06 191.82 NaN 139 184.32 209.17 179.64 204.49 NaN ## get unique ship ave data getUniqueShipsData <- function(df){ results = data.frame() for (thisMMSI in unique(df$mmsi)){ tmp = subset(df, mmsi == thisMMSI & ddB_SL_Emp > 0 & ddB_SL_Emp < 190 & !is.nan(ddB_SL_Emp)) # 190 is an outlier !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! count = nrow(tmp) if (count > 0){ meanEmp = mean(tmp$ddB_SL_Emp) meanSOG = mean(tmp$sog) medianEmp = median(tmp$ddB_SL_Emp) stdevEmp<-sd(tmp$ddB_SL_Emp) meanEmpAbs = mean(tmp$ddB_SL_EmpAbs) medianEmpAbs = median(tmp$ddB_SL_EmpAbs) stdevEmpAbs<-sd(tmp$ddB_SL_EmpAbs) thisClass = tmp$classA[1] results = rbind(results,data.frame(thisMMSI,count,meanEmp,medianEmp,stdevEmp,meanSOG,meanEmpAbs,medianEmpAbs,stdevEmpAbs,thisClass)) } } results } uniqueShips = getUniqueShipsData(data_1_3_BB_100) uniqueShipsSOG = getUniqueShipsData(data_1_3_BB_100) uniqueShipsSorted = uniqueShips[order(uniqueShips$meanEmp),] tail(uniqueShipsSorted) thisMMSI count meanEmp medianEmp stdevEmp meanEmpAbs medianEmpAbs stdevEmpAbs thisClass 1083 311288000 1 185.33 185.33 NA 185.63 185.63 NA Cargo 1072 235768000 1 185.35 185.35 NA 203.03 203.03 NA Cargo 1148 416248000 1 185.71 185.71 NA 196.24 196.24 NA Container ship 1085 372119000 1 185.72 185.72 NA 195.42 195.42 NA Bulk carrier 443 241085000 1 186.20 186.20 NA 186.20 186.20 NA Bulk carrier 1101 241145000 1 186.72 186.72 NA 194.38 194.38 NA Bulk carrier head(uniqueShipsSorted) thisMMSI count meanEmp medianEmp stdevEmp meanEmpAbs medianEmpAbs stdevEmpAbs thisClass 1505 367323180 1 135.26 135.26 NA 135.26 135.26 NA Passenger 357 316017734 1 139.63 139.63 NA 140.48 140.48 NA Miscellaneous 1057 316015742 1 140.65 140.65 NA 141.94 141.94 NA Miscellaneous 1013 303454000 1 142.64 142.64 NA 153.06 153.06 NA Pleasure craft 1419 367473060 1 142.70 142.70 NA 146.19 146.19 NA Pleasure craft 901 367027560 1 144.01 144.01 NA 144.01 144.01 NA Pleasure craft tst=subset(uniqueShipsSorted,count>5) > head(tst) thisMMSI count meanEmp medianEmp stdevEmp meanEmpAbs medianEmpAbs stdevEmpAbs thisClass 570 316013411 9 150.8822 147.020 13.36904 151.0133 147.030 13.42234 Military 113 316011515 8 157.9737 155.155 11.76059 159.3050 155.320 14.15541 Military 247 316013412 6 158.7133 156.345 10.22294 159.7350 156.570 10.08578 Military 207 316192000 6 158.8867 161.405 16.65098 159.0267 161.535 16.67285 Military 73 316010415 16 159.2144 159.420 11.32640 159.8675 162.430 11.52905 Military 211 316009418 8 159.2612 162.300 11.14113 159.4400 162.390 11.11161 Military ######################################### construct power and accumulated power buildPwr_CumPwr <- function(df){ df$pwr = 10^(df$meanEmp/10.0) dfSorted=data.frame() dfSorted = df[order(df$pwr),] N=length(df$pwr) dfSorted$cumPwr = 1 dfSorted$cumDb = 1 dfSorted$index[1]=1; dfSorted$cumPwr[1] = dfSorted$pwr[1] for (i in 2:N){ dfSorted$cumPwr[i] = dfSorted$cumPwr[i-1]+dfSorted$pwr[i] dfSorted$index[i] = i; } dfSorted$cumDb = 10*log10(dfSorted$cumPwr) dfSorted$cumPwrFrac = dfSorted$cumPwr/max(dfSorted$cumPwr) dfSorted } uniqueShipsSortedPwr = buildPwr_CumPwr(uniqueShipsSorted) head(uniqueShipsSortedPwr) thisMMSI count meanEmp medianEmp stdevEmp meanEmpAbs medianEmpAbs stdevEmpAbs thisClass pwr cumPwr cumDb index 1505 367323180 1 135.26 135.26 NA 135.26 135.26 NA Passenger 3.357376e+13 3.357376e+13 135.2600 1 357 316017734 1 139.63 139.63 NA 140.48 140.48 NA Miscellaneous 9.183326e+13 1.254070e+14 140.9832 2 1057 316015742 1 140.65 140.65 NA 141.94 141.94 NA Miscellaneous 1.161449e+14 2.415519e+14 143.8301 3 1013 303454000 1 142.64 142.64 NA 153.06 153.06 NA Pleasure craft 1.836538e+14 4.252057e+14 146.2860 4 1419 367473060 1 142.70 142.70 NA 146.19 146.19 NA Pleasure craft 1.862087e+14 6.114144e+14 147.8634 5 901 367027560 1 144.01 144.01 NA 144.01 144.01 NA Pleasure craft 2.517677e+14 8.631821e+14 149.3610 6 tail(uniqueShipsSortedPwr) 1083 311288000 1 185.33 185.33 NA 185.63 185.63 NA Cargo 3.411929e+18 6.515848e+20 208.1397 1589 1072 235768000 1 185.35 185.35 NA 203.03 203.03 NA Cargo 3.427678e+18 6.550124e+20 208.1625 1590 1148 416248000 1 185.71 185.71 NA 196.24 196.24 NA Container ship 3.723917e+18 6.587364e+20 208.1871 1591 1085 372119000 1 185.72 185.72 NA 195.42 195.42 NA Bulk carrier 3.732502e+18 6.624689e+20 208.2117 1592 443 241085000 1 186.20 186.20 NA 186.20 186.20 NA Bulk carrier 4.168694e+18 6.666376e+20 208.2389 1593 1101 241145000 1 186.72 186.72 NA 194.38 194.38 NA Bulk carrier 4.698941e+18 6.713365e+20 208.2694 1594 ############################################################################## construct strategies # Source Level vs Index plot ggplot() + geom_line(data=uniqueShipsSortedPwr,aes(x=index,y=meanEmp),size=1, color = 'black') + theme_bw()+ # labs(title="Source level vs Ship order") + xlab("ship number")+ ylab(expression(paste("Sourcelevel (dB re 1 ", mu, "Pa"))) + theme(panel.grid.major = element_line(colour="black", size=0.10), panel.grid.minor.y=element_line(color="black",size=0.05), panel.grid.minor.x=element_line(color="black",size=0.05)) # Create plots plot1.x = qplot(y = meanEmp, x = index, xlab="Ship order (red)",color="red", ylab="Source level (dB re 1 microPa @ 1m)", data = uniqueShipsSortedPwr) plot2.x = qplot(y = meanEmp, x = cumPwrFrac, xlab="Cumulative power fraction (black)",data = uniqueShipsSortedPwr) # Run function ggplot_dual_axis(plot1.x, plot2.x, "x") uniqueShipsSortedPwr$shipFraction = uniqueShipsSortedPwr$index/max(uniqueShipsSortedPwr$index) plot1.x = qplot(y = meanEmp, x = shipFraction, xlab="Fraction of ship population (red)",color="red", ylab="Source level (dB re 1 microPa @ 1m)", data = uniqueShipsSortedPwr) plot2.x = qplot(y = meanEmp, x = cumPwrFrac, xlab="Cumulative power fraction (black)",data = uniqueShipsSortedPwr) # Run function ggplot_dual_axis(plot1.x, plot2.x, "x") #################### here is the current data set summary shipclassSummary(uniqueShipsSortedPwr) #function defined below [1] "Total ships = 1581" count countUnique percent meanEmp medianEmp stdevEmp cls 1 734 734 46.047679 173.0891 173.085 0.3920720 Bulk carrier 2 206 206 12.923463 175.0914 175.415 0.4354632 Cargo 3 206 206 12.923463 178.3301 178.635 0.4849865 Container ship 12 111 111 6.963614 175.4346 175.715 0.4848229 Vehicle carrier 10 101 101 6.336261 174.2614 174.360 0.4393965 Tanker 11 85 85 5.332497 169.0829 169.900 0.5026540 Tug 8 35 35 2.195734 157.5417 156.470 0.2840286 Pleasure craft 7 31 31 1.944793 165.8023 167.410 0.4951376 Passenger 4 28 28 1.756587 161.0272 163.290 0.4973475 Fishing 6 21 21 1.317440 161.1324 163.590 0.4629100 Miscellaneous 5 19 19 1.191970 161.3436 161.065 0.3746343 Military 9 4 4 0.250941 167.7547 166.985 0.5773503 Research ############ sink enough ships to get 1/2 total power idxMax = max(uniqueShipsSortedPwr$index) idx3db = min(subset(uniqueShipsSortedPwr,cumPwrFrac>0.5)$index) fracRemoved = (idxMax - idx3db)/idxMax [1] 0.1505646 print(uniqueShipsSortedPwr$meanEmp[idx3db]) [1] 178.9167 shipsControlled_15Percent = subset(uniqueShipsSortedPwr, index>idx3db) shipclassSummary(shipsControlled_15Percent) ########## function defined below count countUnique percent meanEmp medianEmp stdevEmp cls 1 94 94 39.1666667 180.9853 180.365 0.3778347 Bulk carrier 3 90 90 37.5000000 181.2325 181.040 0.4875478 Container ship 2 36 36 15.0000000 181.1119 180.585 0.3779645 Cargo 5 12 12 5.0000000 180.9831 180.815 0.3892495 Tanker 7 5 5 2.0833333 181.2720 180.960 0.5477226 Vehicle carrier 6 2 2 0.8333333 179.2500 179.250 0.0000000 Tug 4 1 1 0.4166667 179.8800 179.880 NA Miscellaneous ############# - 5 dB strategy have loudest ships reduce by 5 dB until total power is 1/2 of original idx = (max(uniqueShipsSortedPwr$index)) cumPwrMax = max(uniqueShipsSortedPwr$cumPwr) cumPwrMod = cumPwrMax while(cumPwrMod > cumPwrMax/2){ cumPwrMod = cumPwrMod - (1 - .3162) * uniqueShipsSortedPwr$pwr[idx] # idea is that each loudest ship removes 31.6% of its power from total idx = idx - 1 } print(idx) [1] 1082 print(uniqueShipsSortedPwr$meanEmp[idx]) [1] 176.33 theFracAffected = (max(uniqueShipsSortedPwr$index) - idx)/max(uniqueShipsSortedPwr$index) print(theFracAffected) [1] 0.3212045 shipsControlled_minus5 = subset(uniqueShipsSortedPwr, index>idx) shipclassSummary(shipsControlled_minus5) ########## function defined below count countUnique percent meanEmp medianEmp stdevEmp cls 1 186 186 36.3281250 179.2705 178.9850 0.4156556 Bulk carrier 3 154 154 30.0781250 179.8000 179.5225 0.4826152 Container ship 2 88 88 17.1875000 178.9538 178.3200 0.4479140 Cargo 8 41 41 8.0078125 177.7484 177.3900 0.4800915 Vehicle carrier 6 34 34 6.6406250 178.6961 178.2050 0.3869530 Tanker 7 5 5 0.9765625 178.4170 178.3900 0.4472136 Tug 4 1 1 0.1953125 179.8800 179.8800 NA Miscellaneous 5 1 1 0.1953125 176.4000 176.4000 NA Passenger ############# - 10 dB strategy have loudest ships reduce by 10 dB until total power is 1/2 of original idx = (max(uniqueShipsSortedPwr$index)) cumPwrMax = max(uniqueShipsSortedPwr$cumPwr) cumPwrMod = cumPwrMax while(cumPwrMod > cumPwrMax/2){ cumPwrMod = cumPwrMod - 0.9 * uniqueShipsSortedPwr$pwr[idx] # idea is that each loudest ship removes 90% of its power from total idx = idx - 1 } print(idx) #[1] 1301 print(uniqueShipsSortedPwr$meanEmp[idx]) [1] 178.3 theFracAffected = (max(uniqueShipsSortedPwr$index) - idx)/max(uniqueShipsSortedPwr$index) print(theFracAffected) [1] 0.1838143 shipsControlled_minus10 = subset(uniqueShipsSortedPwr, index>idx) shipclassSummary(shipsControlled_minus10) ########## function defined below count countUnique percent meanEmp medianEmp stdevEmp cls 3 112 112 38.2252560 180.7268 180.385 0.4785455 Container ship 1 110 110 37.5426621 180.6248 180.160 0.3874598 Bulk carrier 2 44 44 15.0170648 180.6432 180.210 0.3699894 Cargo 5 15 15 5.1194539 180.5051 180.300 0.3518658 Tanker 7 7 7 2.3890785 180.4521 179.490 0.5345225 Vehicle carrier 6 4 4 1.3651877 178.8075 178.775 0.0000000 Tug 4 1 1 0.3412969 179.8800 179.880 NA Miscellaneous ############################################ Move em all down enough to get the -3 dB # find pwr of ship at the half total power point moveShipsDownOne <- function(df){ idx = (max(df$index)) cumPwrMax = max(df$cumPwr) print(sprintf("num=%4.0f",idx)) Nmoved=0 cumPwrMod = cumPwrMax while(cumPwrMod > cumPwrMax/2){ Nmoved = Nmoved + 1 cumPwrMod = cumPwrMod - df$pwr[idx] + # idea is that each loudest ship drops SL to next quieter one Nmoved * df$pwr[idx-1] # and this reduces the total power by the difference print(sprintf("idx=%4.0f totPwr=%e pwr%e Nmoved=%4.0f",idx,cumPwrMod,df$pwr[idx],Nmoved)) df$pwr[idx-1] = df$pwr[idx-1] + df$pwr[idx-1]*Nmoved #move these Nmoved down one and add to pwr idx = idx - 1 } print(idx) } moveShipsDownOne <- function(df){ idxMax = (max(df$index)) cumPwrMax = max(df$cumPwr) print(sprintf("idxMax=%4.0f",idxMax)) Nmoved=0 cumPwrMod = cumPwrMax while(cumPwrMod > cumPwrMax/2){ Nmoved = Nmoved + 1 iStart = idxMax-Nmoved iStop = idxMax-1 for (i in iStart:iStop){ deltaPwr = df$pwr[i+1] - df$pwr[i] df$pwr[i+1] = df$pwr[i] cumPwrMod = cumPwrMod - deltaPwr } # print(sprintf("%4.0f %e",Nmoved,cumPwrMod)) } idxThresh = idxMax-Nmoved percent = 100*(idxMax-Nmoved)/idxMax print(sprintf("Nmoved=%4.0f idxMax=%4.0f index=%4.0f percent=%4.3f dBthreshold=%4.2f ",Nmoved, idxMax,idxThresh,percent,df$meanEmp[idxThresh])) df$modSL = 10*log10(df$pwr) df$cumPwr[1] = df$pwr[1] for (i in 2:idxMax) df$cumPwr[i] = df$cumPwr[i-1] + df$pwr[i] df } tst = moveShipsDownOne(uniqueShipsSortedPwr) [1] "idxMax=1594" [1] "Nmoved= 677 idxMax=1594 index= 917 percent=57.528 dBthreshold=175.05 " > tail(tst) thisMMSI count meanEmp medianEmp stdevEmp meanEmpAbs medianEmpAbs stdevEmpAbs thisClass pwr cumPwr cumDb 1083 311288000 1 185.33 185.33 NA 185.63 185.63 NA Cargo 3.198895e+17 3.338628e+20 208.1397 1072 235768000 1 185.35 185.35 NA 203.03 203.03 NA Cargo 3.198895e+17 3.341827e+20 208.1625 1148 416248000 1 185.71 185.71 NA 196.24 196.24 NA Container ship 3.198895e+17 3.345026e+20 208.1871 1085 372119000 1 185.72 185.72 NA 195.42 195.42 NA Bulk carrier 3.198895e+17 3.348225e+20 208.2117 443 241085000 1 186.20 186.20 NA 186.20 186.20 NA Bulk carrier 3.198895e+17 3.351423e+20 208.2389 1101 241145000 1 186.72 186.72 NA 194.38 194.38 NA Bulk carrier 3.198895e+17 3.354622e+20 208.2694 index cumPwrFrac shipFraction modSL 1083 1589 0.9705785 0.9968632 175.05 1072 1590 0.9756842 0.9974906 175.05 1148 1591 0.9812313 0.9981179 175.05 1085 1592 0.9867911 0.9987453 175.05 443 1593 0.9930006 0.9993726 175.05 1101 1594 1.0000000 1.0000000 175.05 shipDif = uniqueShipsSortedPwr$meanEmp - tst$modSL # this is the amout the various ships have to reduce their SL [1] 8.530000 8.540000 8.615000 8.740000 8.740000 8.740000 8.757273 8.864000 8.875000 8.990000 9.020000 9.045000 9.450000 [14] 9.510000 9.510000 9.760000 9.780000 10.160000 10.260000 10.280000 10.300000 10.660000 10.670000 11.150000 11.670000 dif = as.data.frame(shipDif) dif = subset(dif,shipDif>0.1) summary(dif) shipDif Min. : 0.508 1st Qu.: 1.750 Median : 3.127 Mean : 3.612 3rd Qu.: 5.041 Max. :11.670 line1 = data.frame(x=c(0.85,0.85), y=c(140,179)) plt <- ggplot(data=uniqueShipsSortedPwr, aes(x=shipFraction,y=meanEmp)) + geom_line()+ geom_point(aes(x=.85, y=179),size=3)+ geom_line(data=data.frame(x=c(0.85,0.85),y=c(140,179)),aes(x=x,y=y),size=0.3, linetype="dashed")+ # annotate("text", x = 0.86, y = 140, label = "85%") + geom_text(aes(x=.92,y=177),label="179.0")+ geom_point(aes(x=0.57, y=175), size =3)+ geom_line(data=data.frame(x=c(0.57,0.57),y=c(140,175)),aes(x=x,y=y),size=0.3, linetype = "dotted")+ # annotate("text", x = 0.58, y = 140, label = "57%") + geom_text(aes(x=.64,y=173),label="175.4")+ xlab("Fraction of ship population")+ ylab(expression(paste("Source level (dB re 1 ", mu, "Pa @ 1m)"))) + scale_y_continuous(limits=c(140,190)) + theme_bw()+ ## labs(title="Ship SL's (sorted low to high)")+ theme(panel.grid.major = element_line(colour="black", size=0.00),panel.grid.minor.y=element_line(color="black",size=0.00)) filename = "./sortedSLs" ExportPlot(plt,filename,4,3) #################################################################################### look at lengths vs year of manufacture # dataset with lengths and year data_1_3_BB_100 dataLo179 = subset(data_1_3_BB_100, length>0 & year>0 & !is.nan(ddB_SL_Emp) & ddB_SL_Emp < 179) dataHi179 = subset(data_1_3_BB_100, length>0 & year>0 & !is.nan(ddB_SL_Emp) & ddB_SL_Emp >=179) dataLo175 = subset(data_1_3_BB_100, length>0 & year>0 & !is.nan(ddB_SL_Emp) & ddB_SL_Emp < 175) dataHi175 = subset(data_1_3_BB_100, length>0 & year>0 & !is.nan(ddB_SL_Emp) & ddB_SL_Emp >=175.38) shipclassResultsA(dataHi175) [1] "total count = 1093 total unique = 677" count countUnique meanEmp medianEmp stdevEmp meanEmpAbs medianEmpAbs stdevEmpAbs thisclass 1 359 276 178.7164 178.270 2.4997443 183.7591 179.560 10.4186441 Bulk carrier 2 127 89 178.6065 178.140 2.2032456 186.6132 180.960 11.5196360 Cargo 3 421 184 179.6919 179.400 2.7892389 189.7146 185.140 12.0996005 Container ship 4 1 1 177.9800 177.980 NA 178.0300 178.030 NA Fishing 5 3 3 175.9233 175.950 0.4905439 176.2200 176.130 0.2286919 Passenger 6 56 46 178.5021 177.750 2.2656420 185.7961 179.415 12.4106048 Tanker 7 26 12 177.3565 176.800 1.8652462 191.0431 186.270 13.9790199 Tug 8 100 66 177.4405 177.045 1.6902958 186.6603 179.730 13.0260620 Vehicle carrier shipclassResultsA(dataLo175) ggplot()+geom_point(data=dataLo,aes(x=year,y=length), color="black", size = 1) + geom_point(data=dataHi,aes(x=year,y=length), color="red", size = 1) + xlab("Year of manufacture")+ ylab("Length (m)")+ theme_bw()+ labs(title="Red is ship SL's >= 179 dB")+ theme(panel.grid.major = element_line(colour="black", size=0.00),panel.grid.minor.y=element_line(color="black",size=0.00)) ggplot()+geom_point(data=dataLo,aes(x=year,y=length), color="black", size = 1) + geom_point(data=dataHi,aes(x=year,y=length), color="red", size = 1) + xlab("Year of manufacture")+ ylab("Length (m)")+ theme_bw()+ labs(title="Red is ship SL's >= 179 dB")+ theme(panel.grid.major = element_line(colour="black", size=0.00),panel.grid.minor.y=element_line(color="black",size=0.00)) + facet_wrap(~classA) #################################################### speed limit head(uniqueShipsSOG) thisMMSI count meanEmp medianEmp stdevEmp meanSOG meanEmpAbs medianEmpAbs stdevEmpAbs thisClass 1 316011773 5 171.3980 171.93 2.576494 9.62 171.6340 172.000 2.54877 Military 2 247279500 1 173.8800 173.88 NA 14.00 179.5000 179.500 NA Bulk carrier 3 245855000 1 175.0300 175.03 NA 14.30 175.8000 175.800 NA Cargo 4 316005594 3 157.4067 151.35 12.417223 4.40 158.0967 151.390 13.53008 Tug 5 372359000 1 177.0600 177.06 NA 14.20 191.8200 191.820 NA Bulk carrier 6 636091891 10 178.3690 179.16 2.072585 21.10 185.2210 180.455 11.96978 Container ship uniqueShipsSOGSorted = uniqueShipsSOG[order(uniqueShipsSOG$meanEmp),] uniqueShipsSOGSortedPwr = buildPwr_CumPwr(uniqueShipsSOGSorted) speedLimit = 11.8 # speed limit for -3dB overall Nships = length(uniqueShipsSOGSortedPwr$thisMMSI) Ncontrolled = 0 AveSOGreduction = 0 for (i in 1:length(uniqueShipsSOGSortedPwr$meanSOG)){ uniqueShipsSOGSortedPwr$newSOG[i] = min(speedLimit,uniqueShipsSOGSortedPwr$meanSOG[i]) delta = uniqueShipsSOGSortedPwr$meanSOG[i] - uniqueShipsSOGSortedPwr$newSOG[i] uniqueShipsSOGSortedPwr$SOGreduction[i] = delta uniqueShipsSOGSortedPwr$SL_model[i] = uniqueShipsSOGSortedPwr$meanEmp[i] - delta if (delta > 0){ Ncontrolled = Ncontrolled + 1 AveSOGreduction = AveSOGreduction + delta # print(sprintf("i=%4.0f SL_mod=%4.2f meanEmp=%4.2f delta=%4.2f",i,uniqueShipsSOGSortedPwr$SL_model[i],uniqueShipsSOGSortedPwr$meanEmp[i],delta )) } } AveSOGreduction = AveSOGreduction/Ncontrolled print(sprintf("MODEL speedlimit = %4.2f Fraction of ships controlled = %4.3f Average speed reduction (kts) = %4.2f",speedLimit,Ncontrolled/Nships,AveSOGreduction )) print(summary(subset(uniqueShipsSOGSortedPwr,uniqueShipsSOGSortedPwr$SOGreduction>0)$SOGreduction)) model = build_Model_Pwr_CumPwr(uniqueShipsSOGSortedPwr) tail(model) thisMMSI count meanEmp medianEmp stdevEmp meanSOG meanEmpAbs medianEmpAbs stdevEmpAbs thisClass pwr cumPwr 288 636012547 2 183.665 183.665 3.203194 12.35 200.18 200.18 0.5374012 Bulk carrier 2.325412e+18 6.046867e+20 388 308802000 1 183.410 183.410 NA 9.70 183.45 183.45 NA Cargo 2.192805e+18 5.955566e+20 1279 372193000 1 183.790 183.790 NA 12.00 183.90 183.90 NA Bulk carrier 2.393316e+18 6.094733e+20 322 235082802 1 184.500 184.500 NA 11.80 184.52 184.52 NA Bulk carrier 2.818383e+18 6.296747e+20 1083 311288000 1 185.330 185.330 NA 12.00 185.63 185.63 NA Cargo 3.411929e+18 6.515848e+20 443 241085000 1 186.200 186.200 NA 11.30 186.20 186.20 NA Bulk carrier 4.168694e+18 6.666376e+20 cumDb index newSOG SOGreduction SL_model modPwr modcumPwr modcumDb modindex modcumPwrFrac 288 207.8153 1572 11.8 0.55 183.115 2.048802e+18 3.240061e+20 205.1055 1589 0.9565322 388 207.7492 1568 9.7 0.00 183.410 2.192805e+18 3.261989e+20 205.1348 1590 0.9630058 1279 207.8495 1574 11.8 0.20 183.590 2.285599e+18 3.284845e+20 205.1651 1591 0.9697534 322 207.9912 1582 11.8 0.00 184.500 2.818383e+18 3.313029e+20 205.2023 1592 0.9780738 1083 208.1397 1589 11.8 0.20 185.130 3.258367e+18 3.345613e+20 205.2448 1593 0.9876932 443 208.2389 1593 11.3 0.00 186.200 4.168694e+18 3.387300e+20 205.2985 1594 1.0000000 hist(model$SOGreduction) summary(uniqueShipsSOGSortedPwr$meanEmp) Min. 1st Qu. Median Mean 3rd Qu. Max. 135.3 169.8 174.0 173.0 177.3 186.7 summary(uniqueShipsSOGSortedPwr$SL_model) Min. 1st Qu. Median Mean 3rd Qu. Max. 124.6 156.7 160.5 159.7 163.5 174.9 hist(model$meanSOG) require(stats) reg<-lm(meanEmp ~ meanSOG, data = model) plot(model$meanSOG,model$meanEmp) abline(reg) summary(reg) hist(model$newSOG,breaks=100) speedControlledShips <- subset(uniqueShipsSOGSortedPwr,uniqueShipsSOGSortedPwr$SOGreduction>0) shipclassSummary(speedControlledShips) [1] "Total ships = 1308" count countUnique percent meanEmp medianEmp stdevEmp cls 1 655 655 49.6588324 173.1340 173.1150 0.4025794 Bulk carrier 3 206 206 15.6178923 178.3301 178.6350 0.4849865 Container ship 2 183 183 13.8741471 175.2841 175.4000 0.4468375 Cargo 12 110 110 8.3396513 175.5057 175.7175 0.4857434 Vehicle carrier 10 96 96 7.2782411 174.1880 174.3150 0.4411657 Tanker 7 23 23 1.7437453 168.0740 169.2750 0.4869848 Passenger 8 16 16 1.2130402 161.1128 162.7250 0.2500000 Pleasure craft 5 5 5 0.3790751 167.5709 168.4150 0.4472136 Military 6 5 5 0.3790751 159.5115 156.3500 0.4472136 Miscellaneous 11 5 5 0.3790751 173.1940 172.4200 0.0000000 Tug 4 2 2 0.1516300 172.6980 172.9550 0.7071068 Fishing 9 2 2 0.1516300 171.0500 171.0500 0.0000000 Research printPercents <- function(){ controlled <- shipclassSummary(speedControlledShips) totpop <- shipclassResultsA(data_1_3_BB_100) print("Number in population Number controlled Percentage of population class controlled Class") classes = totpop$thisclass for (aclass in classes){ totCount = subset(totpop,thisclass==aclass)$countUnique controlCount = subset(controlled,cls==aclass)$countUnique percent = 100*controlCount/totCount print(sprintf(" %4.0f %4.0f %4.2f %s",totCount,controlCount,percent,aclass)) } } build_Model_Pwr_CumPwr <- function(df){ df$modPwr = 10^(df$SL_model/10.0) dfSorted=data.frame() dfSorted = df[order(df$modPwr),] N=length(df$modPwr) dfSorted$modcumPwr = 1 dfSorted$modcumDb = 1 dfSorted$modindex[1]=1; dfSorted$modcumPwr[1] = dfSorted$modPwr[1] for (i in 2:N){ dfSorted$modcumPwr[i] = dfSorted$modcumPwr[i-1]+dfSorted$modPwr[i] dfSorted$modindex[i] = i; } dfSorted$modcumDb = 10*log10(dfSorted$modcumPwr) dfSorted$modcumPwrFrac = dfSorted$modcumPwr/max(dfSorted$modcumPwr) print(sprintf("ratio of total population power (model/population) %4.4f",dfSorted$modcumPwr[N]/dfSorted$cumPwr[N])) dfSorted } ############################## helper stuff library(gtable) install.packages("gridExtra") library(grid) library(gridExtra) ggplot_dual_axis = function(plot1, plot2, which.axis = "x") { #### add axis that scales with original (x and/or y) ############ https://gist.github.com/jslefche/e4c0e9f57f0af49fca87 # Update plot with transparent panel plot2 = plot2 + theme(panel.background = element_rect(fill = NA)) + ## VV theme(plot.margin = unit(c(0.7, 1.5, 0.4, 0.4), "cm")) grid.newpage() # Increase right margin if which.axis == "y" if(which.axis == "y") plot1 = plot1 + theme(plot.margin = unit(c(0.7, 1.5, 0.4, 0.4), "cm")) # Extract gtable g1 = ggplot_gtable(ggplot_build(plot1)) g2 = ggplot_gtable(ggplot_build(plot2)) # Overlap the panel of the second plot on that of the first pp = c(subset(g1$layout, name == "panel", se = t:r)) g = gtable_add_grob(g1, g2$grobs[[which(g2$layout$name=="panel")]], pp$t, pp$l, pp$b, pp$l) # Steal axis from second plot and modify axis.lab = ifelse(which.axis == "x", "axis-b", "axis-l") ia = which(g2$layout$name == axis.lab) ga = g2$grobs[[ia]] ax = ga$children[[2]] # Switch position of ticks and labels if(which.axis == "x") ax$heights = rev(ax$heights) else ax$widths = rev(ax$widths) ax$grobs = rev(ax$grobs) if(which.axis == "x") ax$grobs[[2]]$y = ax$grobs[[2]]$y - unit(1, "npc") + unit(0.15, "cm") else ax$grobs[[1]]$x = ax$grobs[[1]]$x - unit(1, "npc") + unit(0.15, "cm") # Modify existing row to be tall enough for axis if(which.axis == "x") g$heights[[2]] = g$heights[g2$layout[ia,]$t] # Add new row or column for axis label if(which.axis == "x") { g = gtable_add_grob(g, ax, 2, 4, 2, 4) g = gtable_add_rows(g, g2$heights[1], 1) g = gtable_add_grob(g, g2$grob[[6]], 2, 4, 2, 4) } else { g = gtable_add_cols(g, g2$widths[g2$layout[ia, ]$l], length(g$widths) - 1) g = gtable_add_grob(g, ax, pp$t, length(g$widths) - 1, pp$b) g = gtable_add_grob(g, g2$grob[[7]], pp$t, length(g$widths), pp$b - 1) } # Draw it grid.draw(g) } ########## demo code # Create fake data.frame data.add.x = data.frame( y1 = runif(100, 0, 100), x1 = runif(100, 0, 100) ) # Add second x-axis that scales with first data.add.x$x2 = (data.add.x$x1 + 50)^0.75 # Create plots plot1.x = qplot(y = y1, x = x1, data = data.add.x) plot2.x = qplot(y = y1, x = x2, data = data.add.x) # Run function ggplot_dual_axis(plot1.x, plot2.x, "x") # Add second y-axis that scales with first data.add.x$y2 = (data.add.x$y^0.5) / 500 # Create plots plot1.y = qplot(y = y1, x = x1, data = data.add.x) plot2.y = qplot(y = y2, x = x1, data = data.add.x) # Run function ggplot_dual_axis(plot1.y, plot2.y, "y") #################################################### get dataframe stats by class -- this is for the unique ships meanEmp files shipclassSummary <- function(df){ results = data.frame() totShips = length(df$thisMMSI); totCount = 0 for (cls in sort(unique(df$thisClass))){ tmp = subset(df, thisClass == cls & meanEmp > 0 & !is.nan(meanEmp)) count = nrow(tmp) if (count > 0){ countUnique = length(unique(tmp$thisMMSI)) totCount = totCount + countUnique percent = 100* countUnique / totShips; meanEmp = mean(tmp$meanEmp) medianEmp = median(tmp$medianEmp) stdevEmp<-sd(!is.na(tmp$stdevEmp)) results = rbind(results,data.frame(count,countUnique,percent, meanEmp,medianEmp,stdevEmp,cls)) } } results = results[order(-results$percent),] print(sprintf("Total ships = %4.0f",totCount)) print(results) } shipclassResultsA <- function(df){ ######### this one is for the ddB_SL_Emmp styled dataframes results = data.frame() total = 0 totalUnique = 0 classes = sort(unique(df$classA)) for (thisclass in classes){ print(thisclass) tmp = subset(df, classA == thisclass & ddB_SL_Emp > 0) count = nrow(tmp) total=total+count countUnique = length(unique(tmp$mmsi)) totalUnique = totalUnique + countUnique meanEmp = mean(tmp$ddB_SL_Emp) medianEmp = median(tmp$ddB_SL_Emp) stdevEmp<-sd(tmp$ddB_SL_Emp) meanEmpAbs = mean(tmp$ddB_SL_EmpAbs) medianEmpAbs = median(tmp$ddB_SL_EmpAbs) stdevEmpAbs<-sd(tmp$ddB_SL_EmpAbs) results = rbind(results,data.frame(count,countUnique,meanEmp,medianEmp,stdevEmp,meanEmpAbs,medianEmpAbs,stdevEmpAbs,thisclass)) } print(sprintf("total transits = %4.0f total unique = %4.0f",total, totalUnique)) results }