![\"\"](\"http:\/\/science.ekqvist.fi\/blogi\/wp-content\/uploads\/2012\/07\/energybalance.jpg\" "\\"energybalance\\"")
<\/a><\/p>\nAs a data source I will use\u00a0 Central Finlands energy balance 2010 data.<\/p>\n
Source of data you will find from this address: Keski-Suomen energiatase 2010<\/a><\/p>\n Ok let’s start to rock with R<\/p>\n To produce Sankey with R we use Function SankeyR (you will find same source as I use in this example from here<\/a>)<\/p>\n We run first SankeyR function. Copy following code into your RGui.<\/p>\n <\/p>\n <\/p>\n After that “Gui-run” we use Central Finlands energy balance data from year 2010…<\/p>\n Result is ok, but just now I have no idea how to produce “sub-branches” of losses. Anyway very nice routine to produce Sankey….And here you will find Sankey diagram in pdf format energybalance<\/a>.<\/p>\n Have fun, Visualize energy balance data As energy engineer I have sometimes need to visualize energy balance data using Sankey diagram. In this example I will introduce how to do that with R. As a data source I will use\u00a0 Central Finlands energy balance 2010 data. Source of data you will find from this address: Keski-Suomen energiatase […] Read More<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"gallery","meta":[],"categories":[19,20,23],"tags":[30,48,51],"_links":{"self":[{"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/posts\/151"}],"collection":[{"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/comments?post=151"}],"version-history":[{"count":0,"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/posts\/151\/revisions"}],"wp:attachment":[{"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/media?parent=151"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/categories?post=151"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/science.ekqvist.fi\/blogi\/wp-json\/wp\/v2\/tags?post=151"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}SankeyR <- function(inputs, losses, unit, labels, format=\"plot\"){\n ########################\n # SankeyR version 1.01 (updated August 10, 2010)\n # is a function for creating Sankey Diagrams in R.\n # See http:\/\/www.sankey-diagrams.com for excellent examples of Sankey Diagrams.\n #\n # OPTIONS:\n # 'inputs' is a vector of input values\n # 'losses' is a vector of loss values\n # 'unit' is a string of the unit\n # 'labels' is a vector of the labels for inputs and losses\n # 'format' is the type of plotting:\n #\u00a0\u00a0 The default is \"plot,\" which produces a plot in the R graphics device.\n #\u00a0\u00a0 Current alternate options include \"pdf\" and \"bmp,\" which produce\n #\u00a0\u00a0 those file types under the name \"Sankey.xxx\" in the current directory.\n #\n # Inputs do not need to equal losses.\u00a0 Any difference will be displayed\n # as a discrepancy in the height of the left and right sides of the diagram.\n # This capability enables the developer to examine imbalances in flows.\n # Percentages are a proportion of the inputs (so, the outputs might not equal 100%).\n #\n # EXAMPLE:\n # Try using these values for the global carbon cycle, from Schlesinger (1997):\n #\u00a0 inputs = c(120,92)\n #\u00a0 losses = c(45,75,90,1,6)\n #\u00a0 unit = \"GtC\/yr\"\n #\u00a0 labels = c(\"GPP\",\"Ocean assimilation\",\"Ra\",\"Rh\",\"Ocean loss\",\"LULCC\",\"Fossil fuel emissions\")\n #\u00a0 SankeyR(inputs,losses,unit,labels)\n #\n # UPDATES:\n # 8\/10\/10 - Added drawing for only one input.\n #\n # CREDITS:\n # Created for R by Aaron BERDANIER\n # send questions or comments to\n # aaron.berdanier@gmail.com\n #\n # SankeyR is based strongly on drawSankey for Matlab,\n # from James SPELLING, KTH-EGI-EKV (spelling@kth.se)\n # http:\/\/leniwiki.epfl.ch\/index.php\/DrawSankey\n #\n # Distributed under Creative Commons Attribution Non-Commercial.\n # Licensees may copy, distribute, display, and perform the work and make\n # derivative works based on it only for noncommercial purposes.\n #\n # Aaron would appreciate notification if you modify or improve this function.\n ########################<\/pre>\n# Calculate fractional losses and inputs\n frLosses = losses\/sum(inputs)\n frInputs = inputs\/sum(inputs)<\/pre>\n
# First input and last output labels\n inputLabel = paste(labels[1],\": \",inputs[1],\" \",unit,\" (\",round(100*frInputs[1],digits=1),\"%)\",sep=\"\")\n lossLabel = paste(labels[length(labels)],\": \",losses[length(losses)],\" \",unit,\" \n(\",round(100*frLosses[length(losses)],digits=1),\"%)\",sep=\"\")<\/pre>\n
########################\n # Calculate position of plot axes (repeat of annotated code below; alternative: save values as vectors...)\n limTop = frInputs[1]; posTop = 0.4; maxy=0\n limBot = 0; posBot = 0.1\n if(length(inputs)>1){\n for(j in 2:length(inputs)){\n rI = max(0.07, abs(frInputs[j]\/2))\n rE = rI + abs(frInputs[j])\n newPosB = posBot + rE*sin(pi\/4) + 0.01\n posBot = newPosB\n arcEx = posBot - rE*sin(seq(0,pi\/4,length.out=100))\n arcEy = limBot - rE*(1-cos(seq(0,pi\/4,length.out=100)))\n arcIx = posBot - rI*sin(seq(0,pi\/4,length.out=100))\n arcIy = limBot - rE + rI*cos(seq(0,pi\/4,length.out=100))\n phiTip = pi\/4 - 2*min(0.05, 0.8*abs(frInputs[j]))\/(rI + rE)\n xTip = posBot - (rE+rI)*sin(phiTip)\/2\n yTip = limBot - rE + (rE+rI)*cos(phiTip)\/2\n limBot = limBot - frInputs[j]\n }\n }else{}\n posTop = posBot + 0.4\n for(i in 1:(length(losses)-1)){\n rI = max(0.07, abs(frLosses[i]\/2))\n rE = rI + abs(frLosses[i])\n arcIx = posTop + rI*sin(seq(0,pi\/2,length.out=100))\n arcIy = limTop + rI*(1 - cos(seq(0,pi\/2,length.out=100)))\n arcEx = posTop + rE*sin(seq(0,pi\/2,length.out=100));\n arcEy = (limTop + rI) - rE*cos(seq(0,pi\/2,length.out=100))\n arEdge = max(0.015, rI\/3)\n arTop\u00a0 = max(0.04, 0.8*frLosses[i])\n arX = posTop + rI + c(0,-arEdge,frLosses[i]\/2,frLosses[i]+arEdge,frLosses[i])\n arY = limTop + rI + c(0,0,arTop,0,0)\n if(max(arY)>maxy){maxy = max(arY)}else{maxy=maxy}\n limTop = limTop - frLosses[i]\n newPos = posTop + rE + 0.01\n posTop = newPos\n }\n newPos = max(posTop, posBot) + max(0.05*limTop,0.05)\n newPos = newPos + 0.8*(limTop-limBot)<\/pre>\nmaxx = newPos\n miny = (limTop-frLosses[length(frLosses)])-max(0.015,abs(frLosses[length(frLosses)]\/4))\n maxy = maxy*2\n minx = 0<\/pre>\n
########################\n # Graphics type?\n if(format!=\"plot\"){\n # Call graphics device\n plottype = switch(format,\n \"pdf\" = pdf(\"Sankey.pdf\", width=11, height=min(8.5,11*(maxy-miny)\/((maxx+3)-minx))),\n \"bmp\" = bmp(\"Sankey.bmp\", width=800*((maxx+3)-minx)\/(maxy-miny), height=800, unit=\"px\",res=144)\n )\n }<\/pre>\n# Create plotting window\n par(mar=c(0,0,0,0),oma=c(0,0,0,0))\n plot(0,0,type=\"n\",xlim=c(-1.5,maxx+1.5),ylim=c(miny,maxy),xaxt=\"n\",yaxt=\"n\")<\/pre>\n
w = 1 # line width<\/pre>\n
# Calculate fractional losses and inputs\n frLosses = losses\/sum(inputs)\n frInputs = inputs\/sum(inputs)<\/pre>\n
# Draw back edge of first input arrow\n lines(c(0.1,0,0.05,0,0.4), c(0,0,frInputs[1]\/2,frInputs[1],frInputs[1]),lwd=w)<\/pre>\n
# First input label\n inputLabel = paste(labels[1],\": \",inputs[1],\" \",unit,\" (\",round(100*frInputs[1],digits=1),\"%)\",sep=\"\")\n fontsize = max(0.5,frInputs[1]*2.5)\n text(0, frInputs[1]\/2, inputLabel, cex=fontsize, pos=2) # try pos=4<\/pre>\n
# Set initial position for the top of the arrows\n limTop = frInputs[1]; posTop = 0.4; maxy=0<\/pre>\n
# set initial position for the bottom of the arrows\n limBot = 0; posBot = 0.1<\/pre>\n
###\n # DRAW ARROWS FOR ADDITIONAL INPUTS\n if(length(inputs)>1){\n for(j in 2:length(inputs)){<\/pre>\n# determine inner and outer arrow radii\n rI = max(0.07, abs(frInputs[j]\/2))\n rE = rI + abs(frInputs[j])\n # push separation point forwards\n newPosB = posBot + rE*sin(pi\/4) + 0.01\n lines(c(posBot,newPosB), c(limBot,limBot), lwd=w)\n posBot = newPosB\n # determine points on the external arc\n arcEx = posBot - rE*sin(seq(0,pi\/4,length.out=100))\n arcEy = limBot - rE*(1-cos(seq(0,pi\/4,length.out=100)))\n # determine points on the internal arc\n arcIx = posBot - rI*sin(seq(0,pi\/4,length.out=100))\n arcIy = limBot - rE + rI*cos(seq(0,pi\/4,length.out=100))\n # draw internal and external arcs\n lines(arcIx, arcIy, lwd=w)\n lines(arcEx, arcEy, lwd=w)<\/pre>\n
# determine arrow point tip\n phiTip = pi\/4 - 2*min(0.05, 0.8*abs(frInputs[j]))\/(rI + rE)\n xTip = posBot - (rE+rI)*sin(phiTip)\/2\n yTip = limBot - rE + (rE+rI)*cos(phiTip)\/2\n # draw back edge of additional input arrows\n lines(c(min(arcEx),xTip,min(arcIx)), c(min(arcEy),yTip,min(arcIy)), lwd=w)<\/pre>\n
# Draw label\n phiText = pi\/2-2*min(0.05,0.8*abs(frInputs[j]))\/(rI+rE)\n xText = posBot-(rE+rI)*sin(phiText)\/3\n yText = limBot-rE\/1.5+(rE+rI)*cos(phiText)\/2\n fullLabel = paste(labels[j],\": \",inputs[j],\" \",unit,\" (\",round(100*frInputs[j],digits=1),\"%)\",sep=\"\")\n fontsize = max(0.5,frInputs[j]*2.5)\n text(xText, yText, fullLabel, cex=fontsize, pos=2)<\/pre>\n
# save new bottom end of arrow\n limBot = limBot - frInputs[j]\n }<\/pre>\n
posTop = posBot + 0.4<\/pre>\n
lines(c(0.4,posTop), c(frInputs[1],frInputs[1]), lwd=w)<\/pre>\n
lines(c(posBot,posBot+(posTop-posBot)\/2), c(limBot,limBot),lwd=w)<\/pre>\n
posMid=posBot+(posTop-posBot)\/2<\/pre>\n
}else{\n lines(c(posBot,posBot+(posTop-posBot)\/2), c(limBot,limBot),lwd=w)\n posMid=posBot+(posTop-posBot)\/2\n }<\/pre>\n###\n # DRAW ARROWS OF LOSSES\n for(i in 1:(length(losses)-1)){<\/pre>\n# Determine inner and outer arrow radii\n rI = max(0.07, abs(frLosses[i]\/2))\n rE = rI + abs(frLosses[i])\n # Determine points on the internal arc\n arcIx = posTop + rI*sin(seq(0,pi\/2,length.out=100))\n arcIy = limTop + rI*(1 - cos(seq(0,pi\/2,length.out=100)))\n # Determine points on the internal arc\n arcEx = posTop + rE*sin(seq(0,pi\/2,length.out=100));\n arcEy = (limTop + rI) - rE*cos(seq(0,pi\/2,length.out=100))\n # Draw internal and external arcs\n lines(arcIx, arcIy, lwd=w)\n lines(arcEx, arcEy, lwd=w)<\/pre>\n
# Determine arrow tip dimensions\n arEdge = max(0.015, rI\/3)\n arTop\u00a0 = max(0.04, 0.8*frLosses[i])\n # Determine points on arrow tip\n arX = posTop + rI + c(0,-arEdge,frLosses[i]\/2,frLosses[i]+arEdge,frLosses[i])\n arY = limTop + rI + c(0,0,arTop,0,0)\n if(max(arY)>maxy){maxy = max(arY)}else{maxy=maxy}\n # Draw tip of losses arrow\n lines(arX, arY, lwd=w)<\/pre>\n# Draw label\n txtX = posTop + rI + frLosses[i]\/2\n txtY = limTop + rI + arTop + 0.05\n fullLabel = paste(labels[i+length(inputs)],\": \",losses[i],\" \",unit,\" \n(\",round(100*frLosses[i],digits=1),\"%)\",sep=\"\")\n fontsize = max(0.5,frLosses[i]*2.5)\n text(txtX, txtY, fullLabel, cex=fontsize, pos=4, srt=35)<\/pre>\n
# Save new position of arrow top\n limTop = limTop - frLosses[i]\n # Advance to new separation point\n newPos = posTop + rE + 0.01\n # Draw top line to new separation point\n lines(c(posTop,newPos), c(limTop,limTop), lwd=w)\n # Save new advancement point\n posTop = newPos<\/pre>\n
# SEPARATION LINES - not implemented yet<\/pre>\n
}<\/pre>\n
###\n # Push the arrow forwards a little after all side-arrows drawn\n newPos = max(posTop, posBot) + max(0.05*limTop,0.05)\n # Draw lines to this new position\n lines(c(posTop,newPos), c(limTop,limTop),lwd=w)\n lines(c(posMid,newPos), c(limTop-frLosses[length(frLosses)],limTop-frLosses[length(frLosses)]),lwd=w)<\/pre>\n
# Draw final arrowhead for the output\n lines(c(newPos,newPos,newPos+max(0.04,0.8*(frLosses[length(frLosses)])),newPos,newPos), \nc(limTop,limTop+max(0.015,abs(frLosses[length<\/pre>\n
(frLosses)]\/6)),limTop-frLosses[length(frLosses)]\/2,(limTop-frLosses[length(frLosses)])-max(0.015,\nabs(frLosses[length(frLosses)]\/6)),(limTop-frLosses<\/pre>\n
[length(frLosses)])),\u00a0 lwd=w)<\/pre>\n
# Save final tip position\n newPos = newPos + 0.8*(frLosses[length(frLosses)])<\/pre>\n
# Last loss label\n lossLabel = paste(labels[length(labels)],\": \",losses[length(losses)],\" \",unit,\" \n(\",round(100*frLosses[length(losses)],digits=1),\"%)\",sep=\"\")\n fontsize = max(0.5,frLosses[length(losses)]*2.5)\n text(newPos+0.05, limTop-frLosses[length(frLosses)]\/2, lossLabel, cex=fontsize, pos=4) # try pos=4<\/pre>\n
# Draw mid-line\n if(limBot<(limTop-frLosses[length(frLosses)])){\n lines(c(posMid,posMid), c(frInputs[1],limBot),lty=2)\n }else{\n lines(c(posMid,posMid), c(frInputs[1],limTop-frLosses[length(frLosses)]),lty=2)\n }<\/pre>\nif(format!=\"plot\"){\n # Close graphics device\n dev.off()\n }}<\/pre>\n#Central Finlands energy balance 2010 (Keski-Suomen energiatase 2010 Source: http:\/\/www.kesto.fi)\n #\n # inputs\n #Oil 4.2 TWh\n #Coal 0,04\n #Peat 3.2\n #Wood fuel 4.1\n #Black ? 2.5\n #REF + other 0.12\n #Water 0.13\n #Elec.import 4.3<\/pre>\n
# lossess\n #Industry 9.5\n #--> Elec 4,65\n #--> process heat 4,85<\/pre>\n
#Buildings 5.1\n #--> district heating 2.55\n #--> wood 0.663\n #--> oil 0.969\n #--> elec 0,867\n #Other elec consumption 1.5\n #Traffic 2.6\n #--> Bensin\u00a0 1.196\n #--> Diesel 1.404\n #total 18.6 TWh<\/pre>\n
inputs = c(4.2,0.04,3.2,4.1,2.5,0.12,4.3,0.13)\n losses = c(9.5,5.1,1.5,2.6)\n unit = \"TWh\/yr\"\n labels = c(\"Oil\",\"Coal\",\"Peat\",\"Wood\",\"Black liquor\",\"REF + other\",\"Elec.import\",\"Water\", \n\"Industry\",\"Buildings\", \"Other elec. consumption\", \"Traffic\")\n SankeyR(inputs,losses,unit,labels)<\/pre>\n
#saving into file jpg. note: picture quality not so good\n jpeg(\"...\/energybalance.jpg\")\n SankeyR(inputs,losses,unit,labels)\n dev.off()<\/pre>\n
#saving into file note: picture quality not so good\n png(\"...\/energybalance.png\")\n SankeyR(inputs,losses,unit,labels)\n dev.off()<\/pre>\n
#saving into file note: note: picture quality not so good\n bmp(\"...\/energybalance.bmp\")\n SankeyR(inputs,losses,unit,labels)\n dev.off()<\/pre>\n
#saving into file note: picture quality excellent\n pdf(\"...\/energybalance.pdf\")\n SankeyR(inputs,losses,unit,labels)\n dev.off()<\/pre>\n
\nMarko<\/p>\n","protected":false},"excerpt":{"rendered":"