# R.E. Benestad # This is a heuristic description of windspeed assuming a Weibull distribution. # We assume some values for the scale and shape parameters: scale=6 shape=2 # Here h is a sample histogram h <-list(breaks=0:20, counts=c(2257,995,4737,3369,4301,2358,2228,844,453,440,142,45,43,6,3,4,1,0,0,0), density=c(1.015477e-01, 4.476739e-02, 2.131288e-01, 1.515792e-01, 1.935121e-01, 1.060920e-01, 1.002430e-01, 3.797354e-02, 2.038154e-02, 1.979663e-02, 6.388914e-03, 2.024656e-03, 1.934671e-03, 2.699541e-04, 1.349771e-04, 1.799694e-04, 4.499235e-05, 0.000000e+00, 0.000000e+00, 0.000000e+00), mids=seq(0.5,19.5,by=1)) # Here, x is the windspeed x <- seq(0,90,by=0.3) X <- (1 + (x-1))^2 dx <- diff(x)[1] # Case # 1 x11() par(col.axis="white",bty="n",cex=1.5,cex.lab=1.2) y <- dweibull(x/3,scale=scale,shape=shape) y <- y/sum(y*dx) #shape2 <- 2.7; scale2=10.5 #y2 <- dweibull(x/3,scale=scale2,shape=shape2)*X; y2 <- y2/sum(y2*dx) y2 <- y*X; y2 <- y2/sum(y2*dx) plot(h\$mids*3,h\$density,type="n",lwd=1,lty=2,xlim=c(0,90),ylim=c(0,0.07), main="Windspeed fraction associated with TCs",xlab="Wind speed (m/s)",ylab="Frequency") lines(x,y,col="red",lwd=3) polygon(c(x[x>33],33),c(y[x>33],0),col="red",border="red") lines(rep(33,2),c(0,1),lty=2,col="grey") lines(rep(43,2),c(0,1),lty=3,col="grey") lines(rep(50,2),c(0,1),lty=3,col="grey") lines(rep(59,2),c(0,1),lty=3,col="grey") lines(rep(70,2),c(0,1),lty=3,col="grey") text(35,0.002,"A",cex=1.5,font=2) text(15,0.002,"a",cex=1.5,font=2) #lines(x,y2,lwd=3,col="grey") a1 <- round(sum(y[x<33]*dx),3) a2 <- round(sum(y2[x<33]*dx),3) text(5,0.066,paste("a=",a1),col="red",font=2,pos=4) #text(5,0.061,paste("a=",a2),col="grey",font=2,pos=4) dev2bitmap(file="windhist1.png",res=150) # Case # 2 x11() par(col.axis="white",bty="n",cex=1.5,cex.lab=1.2) shape2 <- 1.7; scale2=6.5 plot(h\$mids*3,h\$density,type="n",lwd=1,lty=2,xlim=c(0,90),ylim=c(0,0.07), main="Fewer TCs & higher windspeed",xlab="Wind speed (m/s)",ylab="Frequency") lines(x,y,col="red",lwd=3) polygon(c(x[x>33],33),c(y[x>33],0),col="red",border="red") lines(rep(33,2),c(0,1),lty=2,col="grey") lines(rep(43,2),c(0,1),lty=3,col="grey") lines(rep(50,2),c(0,1),lty=3,col="grey") lines(rep(59,2),c(0,1),lty=3,col="grey") lines(rep(70,2),c(0,1),lty=3,col="grey") text(70,0.002,"A",cex=1.5,font=2) text(15,0.002,"a",cex=1.5,font=2) # Monte-Carlo simulation to see if it's possible to get fewer # TCs defined as the area a1 & a2 and more intensive winds # as described by A1 & A2 n <- 10000 scales <- runif(n,min=floor(scale-2),max=ceiling(scale+2)) shapes <- runif(n,min=floor(shape-2),max=ceiling(shape+2)) a.check <- rep(NA,n) ratio0 <- 1 for (i in 1:n) { Y <- dweibull(x/3,scale=scales[i],shape=shapes[i]); Y <- Y/sum(Y*dx) a1 <- sum(y[x<33]*dx); a2 <- sum(Y[x<33]*dx) A1 <- sum(y[x>50]*dx); A2 <- sum(Y[x>50]*dx) a.check[i] <- (a2 > a1) & (A2 > A1) if ( (a2 > a1) & (A2 > A1) ) { ratio <- (a2 - a1)/(A2 - A1) if (ratio > ratio0) { y3 <- Y; A3 <- A2 ratio0 <- ratio } } } print(table(a.check)) lines(x,y3,lwd=3,col="grey") a1 <- round(sum(y[x<33]*dx),3) a2 <- round(sum(y3[x<33]*dx),3) text(5,0.066,paste("a=",a1),col="red",font=2,pos=4) text(5,0.061,paste("a=",a2),col="grey",font=2,pos=4) text(50,0.066,paste("A=",round(A1*1000,3),"e-3"),col="red",font=2,pos=4) text(50,0.061,paste("A=",round(A3*1000,3),"e-3"),col="grey",font=2,pos=4) dev2bitmap(file="windhist2.png",res=150)