docs/tutorials/18 - AgriPV - Coffee Plantation with Tree Modeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4997ec7a",
   "metadata": {},
   "source": [
    "# 12d - AgriPV: Designing for adecuate crop shading "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393445a3",
   "metadata": {},
   "source": [
    "This journal supports the process of designing a solar panel configuration to appropriately represent ideal shading conditions for coffee production underneath elevated solar panels. \n",
    "\n",
    "The coffee trees would be under and/or in between elevated solar panels (panels would be elevated 6, 8, or 10 ft tall). The light/shade analysis helps determine appropriate panel heights and spacings t0 achieve appropriate shading.  The desired level of shading is maximum of 30% (i.e., 70% of normal, unshaded light). \n",
    "\n",
    "Details:\n",
    "*  The coffee plants are expected to be \\~5 ft tall. (5-6 ft tall and 3 ft wide (<a href=\"https://realgoodcoffeeco.com/blogs/realgoodblog/how-to-grow-a-coffee-plant-at-home#:~:text=However%2C%20you%20must%20keep%20in,tall%20and%203%20feet%20wide\">Reference</a>)\n",
    "*\tLocation: 18.202142, -66.759187; (18°12'07.7\"N 66°45'33.1\"W)\n",
    "*\tDesired area of initial analysis: 400-600 ft2 (37-55 m2)\n",
    "*\tRacking: Fixed-tilt panels\n",
    "*\tPanel size: 3.3 feet x 5.4 feet                                    (1m x 1.64m)\n",
    "*\tAnalysis variations:\n",
    "<ul> <li> a.\tPanel height: would like to examine heights of 6 ft, 8 ft, and 10 ft hub height. \n",
    "<li> b.\tPanel spacing (N/W): would like to look at multiple distances (e.g., 2 ft, 3 ft, 4 ft) </li> \n",
    "<li> c.\tInter-Row spacing (E/W): would like to look at multiple distances (e.g., 2 ft, 3 ft, 4 ft)! </li> \n",
    "\n",
    "\n",
    "Steps on this Journal:\n",
    "<ol>\n",
    "    <li> <a href='#step1'> <u><b>Loop to Raytrace and sample irradiance at where Three would be located </u></b></li>\n",
    "    <li> <a href='#step2'> Calculate GHI for Comparisons </li>\n",
    "        <ul><li> <a href='#step2a'> Option 1: Raytrace of Empty Field  </li></ul>\n",
    "        <ul><li> <a href='#step2b'> Option 2: Weather File </li></ul>\n",
    "    <li> <a href='#step3'> Compile Results </li>\n",
    "    <li> <a href='#step4'> Plot Results</li>\n",
    "    <li> <a href='#step5'> <u><b> Raytrace with Tree Geometry <u></b></li>\n",
    "        <ul><li> <a href='#step5a'>Tree Parameters</li></ul>\n",
    "        <ul><li> <a href='#step5b'>Loop to Raytrace and Sample Irradiance at Each side of the Tree (N, S, E, W)</li></ul>\n",
    "        <ul><li> <a href='#step5c'>Single simulation until MakeOct for Getting a PRETTY IMAGE </li></ul>\n",
    "    <li> <a href='#step6'> Compile Results</li>\n",
    "    <li> <a href='#step7'>  Plot </li>\n",
    "\n",
    "</ol>\n",
    "        \n",
    "\n",
    "        \n",
    "![AgriPV Coffee Trees Simulation](../images_wiki/AdvancedJournals/AgriPV_CoffeeTrees.PNG)\n",
    "\n",
    " \n",
    "While we have HPC scripts to do the below simulation, this journals runs all of the above so it might take some time, as there are 109 combinations of parameters explored"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "064f0e47",
   "metadata": {},
   "outputs": [],
   "source": [
    "import bifacial_radiance\n",
    "import os\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd97e1fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "testfolder = str(Path().resolve().parent.parent / 'bifacial_radiance' / 'TEMP' /  'Tutorial_18')\n",
    "if not os.path.exists(testfolder):\n",
    "    os.makedirs(testfolder)\n",
    "    \n",
    "resultsfolder = os.path.join(testfolder, 'results')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "160eec85",
   "metadata": {},
   "source": [
    "### General Parameters and Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd3730bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "lat = 18.202142\n",
    "lon = -66.759187\n",
    "\n",
    "albedo = 0.25 # Grass value from Torres Molina, \"Measuring UHI in Puerto Rico\" 18th LACCEI \n",
    "              # International Multi-Conference for Engineering, Education, and Technology\n",
    "\n",
    "ft2m = 0.3048\n",
    "\n",
    "# Loops\n",
    "clearance_heights = np.array([6.0, 8.0, 10.0])* ft2m\n",
    "xgaps = np.array([2, 3, 4]) * ft2m\n",
    "Ds = np.array([2, 3, 4]) * ft2m    # D is a variable that represents the spacing between rows, not-considering the collector areas.\n",
    "tilts = [round(lat), 10]\n",
    "\n",
    "x = 1.64\n",
    "y = 1        \n",
    "azimuth = 180\n",
    "nMods = 20\n",
    "nRows = 7\n",
    "numpanels = 1\n",
    "moduletype = 'test-module'\n",
    "hpc = False\n",
    "sim_general_name = 'tutorial_18'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5db1643e",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(os.path.join(testfolder, 'EPWs')):\n",
    "    demo = bifacial_radiance.RadianceObj('test',testfolder)  \n",
    "    epwfile = demo.getEPW(lat,lon)    \n",
    "else:\n",
    "    epwfile = r'EPWs\\PRI_Mercedita.AP.785203_TMY3.epw' "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "168ea11a",
   "metadata": {},
   "source": [
    "<a id='step1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48ab5feb",
   "metadata": {},
   "source": [
    "## 1. Loop to Raytrace and sample irradiance at where Three would be located"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43fb3f84",
   "metadata": {},
   "outputs": [],
   "source": [
    "demo = bifacial_radiance.RadianceObj(sim_general_name,str(testfolder))  \n",
    "demo.setGround(albedo)\n",
    "demo.readWeatherFile(epwfile)\n",
    "demo.genCumSky()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8b17974",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ch in range (0, len(clearance_heights)):\n",
    "    \n",
    "    clearance_height = clearance_heights[ch]\n",
    "    for xx in range (0, len(xgaps)):\n",
    "        \n",
    "        xgap = xgaps[xx]\n",
    "\n",
    "        for tt in range (0, len(tilts)):\n",
    "        \n",
    "            tilt = tilts[tt]\n",
    "            for dd in range (0, len(Ds)):\n",
    "                pitch = y * np.cos(np.radians(tilt))+Ds[dd]\n",
    "\n",
    "                sim_name = (sim_general_name+'_ch_'+str(round(clearance_height,1))+\n",
    "                                '_xgap_'+str(round(xgap,1))+\\\n",
    "                                '_tilt_'+str(round(tilt,1))+\n",
    "                                '_pitch_'+str(round(pitch,1)))\n",
    "\n",
    "                # Coffe plant location at:\n",
    "                coffeeplant_x = (x+xgap)/2\n",
    "                coffeeplant_y = pitch/2\n",
    "\n",
    "                demo.makeModule(name=moduletype, x=x, y=y, xgap = xgap)\n",
    "                sceneDict = {'tilt':tilt,'pitch':pitch,'clearance_height':clearance_height,'azimuth':azimuth, 'nMods': nMods, 'nRows': nRows} \n",
    "                scene = demo.makeScene(moduletype=moduletype,sceneDict=sceneDict, hpc=hpc, radname = sim_name)\n",
    "                octfile = demo.makeOct(octname = demo.basename , hpc=hpc)  \n",
    "                analysis = bifacial_radiance.AnalysisObj(octfile=octfile, name=sim_name)\n",
    "\n",
    "                # Modify sensor position to coffee plant location\n",
    "                frontscan, backscan = analysis.moduleAnalysis(scene=scene, sensorsy=1)\n",
    "                groundscan = frontscan.copy() \n",
    "                groundscan['xstart'] = coffeeplant_x\n",
    "                groundscan['ystart'] = coffeeplant_y\n",
    "                groundscan['zstart'] = 0.05\n",
    "                groundscan['orient'] = '0 0 -1'\n",
    "                analysis.analysis(octfile, name=sim_name+'_Front&Back', frontscan=frontscan, backscan=backscan)\n",
    "                analysis.analysis(octfile, name=sim_name+'_Ground&Back', frontscan=groundscan, backscan=backscan)\n",
    "\n",
    "                "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c433ce0a",
   "metadata": {},
   "source": [
    "<a id='step2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff5f9bf7",
   "metadata": {},
   "source": [
    "## 2. Calculate GHI for Comparisons"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3f3e8b2",
   "metadata": {},
   "source": [
    "<a id='step2a'></a>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9f3f7b6",
   "metadata": {},
   "source": [
    "### Option 1: Raytrace of Empty Field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c3afc43",
   "metadata": {},
   "outputs": [],
   "source": [
    "sim_name = 'EMPTY'\n",
    "demo.makeModule(name=moduletype, x=0.001, y=0.001, xgap = 0)\n",
    "sceneDict = {'tilt':0,'pitch':2,'clearance_height':0.005,'azimuth':180, 'nMods': 1, 'nRows': 1} \n",
    "scene = demo.makeScene(moduletype=moduletype,sceneDict=sceneDict, hpc=hpc, radname = sim_name)\n",
    "octfile = demo.makeOct(octname = demo.basename , hpc=hpc)  \n",
    "analysis = bifacial_radiance.AnalysisObj(octfile=octfile, name=sim_name)\n",
    "frontscan, backscan = analysis.moduleAnalysis(scene=scene, sensorsy=1)\n",
    "emptyscan = frontscan.copy() \n",
    "emptyscan['xstart'] = 3\n",
    "emptyscan['ystart'] = 3\n",
    "emptyscan['zstart'] = 0.05\n",
    "emptyscan['orient'] = '0 0 -1'\n",
    "emptybackscan = emptyscan.copy()\n",
    "emptybackscan['orient'] = '0 0 1'\n",
    "analysis.analysis(octfile, name='_EMPTYSCAN', frontscan=emptyscan, backscan=emptybackscan)\n",
    "\n",
    "resname = os.path.join(resultsfolder, 'irr__EMPTYSCAN.csv')\n",
    "data = pd.read_csv(resname)\n",
    "puerto_rico_Year = data['Wm2Front'][0]\n",
    "print(\"YEARLY TOTAL Wh/m2:\", puerto_rico_Year)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d722641",
   "metadata": {},
   "source": [
    "<a id='step2b'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99947d38",
   "metadata": {},
   "source": [
    "### Option 2: Weather File"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "566360fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indexes for start of each month of interest in TMY3 8760 hours file\n",
    "#starts = [2881, 3626, 4346, 5090, 5835]\n",
    "#ends = [3621, 4341, 5085, 5829, 6550]\n",
    "\n",
    "starts = [metdata.datetime.index(pd.to_datetime('2021-05-01 6:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-06-01 6:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-07-01 6:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-08-01 6:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-09-01 6:0:0 -7'))]\n",
    "\n",
    "ends = [metdata.datetime.index(pd.to_datetime('2021-05-31 18:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-06-30 18:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-07-31 18:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-08-31 18:0:0 -7')),\n",
    "          metdata.datetime.index(pd.to_datetime('2021-09-30 18:0:0 -7'))]\n",
    "\n",
    "ghi_PR=[]\n",
    "for ii in range(0, len(starts)):\n",
    "    start = starts[ii]\n",
    "    end = ends[ii]\n",
    "    ghi_PR.append(demo.metdata.ghi[start:end].sum())\n",
    "puerto_Rico_Monthly = ghi_PR     # Wh/m2\n",
    "puerto_Rico_YEAR = demo.metdata.ghi.sum()  # Wh/m2\n",
    "\n",
    "print(\"Monthly Values May-Sept:\", puerto_Rico_Monthly, \"Wh/m2\")\n",
    "print(\"Year Values\", puerto_Rico_YEAR, \"Wh/m2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af748a77",
   "metadata": {},
   "source": [
    "<a id='step3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21682a45",
   "metadata": {},
   "source": [
    "## 3. Compile Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51dad298",
   "metadata": {},
   "outputs": [],
   "source": [
    "ch_all = []\n",
    "xgap_all = []\n",
    "tilt_all = []\n",
    "pitch_all = []\n",
    "FrontIrrad = []\n",
    "RearIrrad = []\n",
    "GroundIrrad = []\n",
    "\n",
    "for ch in range (0, len(clearance_heights)):\n",
    "    \n",
    "    clearance_height = clearance_heights[ch]\n",
    "    for xx in range (0, len(xgaps)):\n",
    "        \n",
    "        xgap = xgaps[xx]\n",
    "\n",
    "        for tt in range (0, len(tilts)):\n",
    "        \n",
    "            tilt = tilts[tt]\n",
    "            for dd in range (0, len(Ds)):\n",
    "                pitch = y * np.cos(np.radians(tilt))+Ds[dd]\n",
    "\n",
    "                # irr_Coffee_ch_1.8_xgap_0.6_tilt_18_pitch_1.6_Front&Back.csv\n",
    "                sim_name = ('irr_Coffee'+'_ch_'+str(round(clearance_height,1))+\n",
    "                                '_xgap_'+str(round(xgap,1))+\\\n",
    "                                '_tilt_'+str(round(tilt,1))+\n",
    "                                '_pitch_'+str(round(pitch,1))+'_Front&Back.csv')\n",
    "\n",
    "                sim_name2 = ('irr_Coffee'+'_ch_'+str(round(clearance_height,1))+\n",
    "                                '_xgap_'+str(round(xgap,1))+\\\n",
    "                                '_tilt_'+str(round(tilt,1))+\n",
    "                                '_pitch_'+str(round(pitch,1))+'_Ground&Back.csv')\n",
    "\n",
    "                ch_all.append(clearance_height)\n",
    "                xgap_all.append(xgap)\n",
    "                tilt_all.append(tilt)\n",
    "                pitch_all.append(pitch)\n",
    "                data = pd.read_csv(os.path.join(resultsfolder, sim_name))\n",
    "                FrontIrrad.append(data['Wm2Front'].item())\n",
    "                RearIrrad.append(data['Wm2Back'].item())\n",
    "                data = pd.read_csv(os.path.join(resultsfolder, sim_name2))\n",
    "                GroundIrrad.append(data['Wm2Front'].item())\n",
    "\n",
    "ch_all = pd.Series(ch_all, name='clearance_height')\n",
    "xgap_all = pd.Series(xgap_all, name='xgap')\n",
    "tilt_all = pd.Series(tilt_all, name='tilt')\n",
    "pitch_all = pd.Series(pitch_all, name='pitch')\n",
    "FrontIrrad = pd.Series(FrontIrrad, name='FrontIrrad')\n",
    "RearIrrad = pd.Series(RearIrrad, name='RearIrrad')\n",
    "GroundIrrad = pd.Series(GroundIrrad, name='GroundIrrad')\n",
    "\n",
    "df = pd.concat([ch_all, xgap_all, tilt_all, pitch_all, FrontIrrad, RearIrrad, GroundIrrad], axis=1)\n",
    "df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5360bdc",
   "metadata": {},
   "source": [
    "#### Let's calculate some relevant metrics for irradiance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "521f9e24",
   "metadata": {},
   "outputs": [],
   "source": [
    "df[['GroundIrrad_percent_GHI']] = df[['GroundIrrad']]*100/puerto_Rico_YEAR\n",
    "df['FrontIrrad_percent_GHI'] = df['FrontIrrad']*100/puerto_Rico_YEAR\n",
    "df['RearIrrad_percent_GHI'] = df['RearIrrad']*100/puerto_Rico_YEAR\n",
    "df['BifacialGain'] = df['RearIrrad']*0.65*100/df['FrontIrrad']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e42529b",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df['GroundIrrad_percent_GHI'].min())\n",
    "print(df['GroundIrrad_percent_GHI'].max())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "229c1d88",
   "metadata": {},
   "source": [
    "<a id='step4'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21972587",
   "metadata": {},
   "source": [
    "## 4. Plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ae01c65",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0a711fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "tilts_l = list(df['tilt'].unique())\n",
    "ch_l = list(df['clearance_height'].unique())\n",
    "print(tilts_l)\n",
    "print(ch_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbb661ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "for tilt in tilts_l:\n",
    "    for clearance_height in ch_l:\n",
    "        df2=df.loc[df['tilt']==tilts[1]]\n",
    "        df3 = df2.loc[df2['clearance_height']==clearance_heights[2]]\n",
    "        df3['pitch']=df3['pitch'].round(1)\n",
    "        df3['xgap']=df3['xgap'].round(1)\n",
    "\n",
    "        sns.set(font_scale=2) \n",
    "        table = df3.pivot('pitch', 'xgap', 'GroundIrrad_percent_GHI')\n",
    "        ax = sns.heatmap(table, cmap='hot', vmin = 50, vmax= 100, annot=True)\n",
    "        ax.invert_yaxis()\n",
    "        figtitle = 'Clearance Height ' + str(clearance_height/ft2m)+'  ft, Tilt ' + str(tilt) + '$^\\circ$'\n",
    "        plt.title(figtitle)\n",
    "        print(table)\n",
    "        plt.show()\n",
    "        \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a221a6a",
   "metadata": {},
   "source": [
    "<a id='step5'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b75d7365",
   "metadata": {},
   "source": [
    "# 5. Raytrace with Tree Geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b277e449",
   "metadata": {},
   "source": [
    "<a id='step5a'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e5149d1",
   "metadata": {},
   "source": [
    "#### Tree parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4476d3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_albedo = 0.165 # Wikipedia [0.15-0.18]\n",
    "trunk_x = 0.8 * ft2m\n",
    "trunk_y = trunk_x\n",
    "trunk_z = 1 * ft2m\n",
    "\n",
    "tree_x = 3 * ft2m\n",
    "tree_y = tree_x\n",
    "tree_z = 4 * ft2m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abb4f71f",
   "metadata": {},
   "source": [
    "<a id='step5b'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d33965eb",
   "metadata": {},
   "source": [
    "#### Loop to Raytrace and Sample Irradiance at Each side of the Tree (N, S, E, W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98857c7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ch in range (0, len(clearance_heights)):\n",
    "    \n",
    "    clearance_height = clearance_heights[ch]\n",
    "    for xx in range (0, len(xgaps)):\n",
    "        \n",
    "        xgap = xgaps[xx]\n",
    "\n",
    "        for tt in range (0, len(tilts)):\n",
    "        \n",
    "            tilt = tilts[tt]\n",
    "            for dd in range (0, len(Ds)):\n",
    "                pitch = y * np.cos(np.radians(tilt))+Ds[dd]\n",
    "\n",
    "                sim_name = (sim_general_name+'_ch_'+str(round(clearance_height,1))+\n",
    "                                '_xgap_'+str(round(xgap,1))+\\\n",
    "                                '_tilt_'+str(round(tilt,1))+\n",
    "                                '_pitch_'+str(round(pitch,1)))\n",
    "\n",
    "                coffeeplant_x = (x+xgap)/2\n",
    "                coffeeplant_y = pitch\n",
    "\n",
    "                demo.makeModule(name=moduletype, x=x, y=y, xgap = xgap)\n",
    "                sceneDict = {'tilt':tilt,'pitch':pitch,'clearance_height':clearance_height,'azimuth':azimuth, 'nMods': nMods, 'nRows': nRows} \n",
    "                scene = demo.makeScene(moduletype=moduletype,sceneDict=sceneDict, hpc=hpc, radname = sim_name)\n",
    "\n",
    "                # Appending the Trees here\n",
    "                text = ''\n",
    "                for ii in range(0,3):\n",
    "                    coffeeplant_x = (x+xgap)/2 + (x+xgap)*ii\n",
    "                    for jj in range(0,3):\n",
    "                        coffeeplant_y = pitch/2 + pitch*jj\n",
    "                        name = 'tree'+str(ii)+str(jj)\n",
    "                        text += '\\r\\n! genrev Metal_Grey tube{}tree t*{} {} 32 | xform -t {} {} {}'.format('head'+str(ii)+str(jj),tree_z, tree_x/2.0, \n",
    "                                                                                                         -trunk_x/2.0 + coffeeplant_x, \n",
    "                                                                                                           -trunk_x/2.0 + coffeeplant_y, trunk_z)\n",
    "                        text += '\\r\\n! genrev Metal_Grey tube{}tree t*{} {} 32 | xform -t {} {} 0'.format('trunk'+str(ii)+str(jj),trunk_z, trunk_x/2.0, \n",
    "                                                                                                         -trunk_x/2.0 + coffeeplant_x, \n",
    "                                                                                                          -trunk_x/2.0 + coffeeplant_y)\n",
    "              \n",
    "                        customObject = demo.makeCustomObject(name,text)\n",
    "                        demo.appendtoScene(radfile=scene.radfiles, customObject=customObject, text=\"!xform -rz 0\")\n",
    "\n",
    "                octfile = demo.makeOct(octname = demo.basename , hpc=hpc)  \n",
    "                analysis = bifacial_radiance.AnalysisObj(octfile=octfile, name=sim_name)\n",
    "\n",
    "\n",
    "                ii = 1\n",
    "                jj = 1\n",
    "                coffeeplant_x = (x+xgap)/2 + (x+xgap)*ii \n",
    "                coffeeplant_y = pitch/2 + pitch*jj\n",
    "                frontscan, backscan = analysis.moduleAnalysis(scene=scene, sensorsy=1)\n",
    "\n",
    "                treescan_south = frontscan.copy()\n",
    "                treescan_north = frontscan.copy()\n",
    "                treescan_east = frontscan.copy()\n",
    "                treescan_west = frontscan.copy()\n",
    "                \n",
    "                treescan_south['xstart'] = coffeeplant_x\n",
    "                treescan_south['ystart'] = coffeeplant_y  - tree_x/2.0 - 0.05\n",
    "                treescan_south['zstart'] = tree_z\n",
    "                treescan_south['orient'] = '0 1 0'\n",
    "\n",
    "                treescan_north['xstart'] = coffeeplant_x\n",
    "                treescan_north['ystart'] = coffeeplant_y  + tree_x/2.0 + 0.05\n",
    "                treescan_north['zstart'] = tree_z\n",
    "                treescan_north['orient'] = '0 -1 0'\n",
    "\n",
    "                treescan_east['xstart'] = coffeeplant_x + tree_x/2.0 + 0.05\n",
    "                treescan_east['ystart'] = coffeeplant_y \n",
    "                treescan_east['zstart'] = tree_z\n",
    "                treescan_east['orient'] = '-1 0 0'\n",
    "\n",
    "                treescan_west['xstart'] = coffeeplant_x - tree_x/2.0 - 0.05\n",
    "                treescan_west['ystart'] = coffeeplant_y \n",
    "                treescan_west['zstart'] = tree_z\n",
    "                treescan_west['orient'] = '1 0 0'\n",
    "    \n",
    "                groundscan = frontscan.copy() \n",
    "                groundscan['xstart'] = coffeeplant_x\n",
    "                groundscan['ystart'] = coffeeplant_y\n",
    "                groundscan['zstart'] = 0.05\n",
    "                groundscan['orient'] = '0 0 -1'\n",
    "                analysis.analysis(octfile, name=sim_name+'_North&South', frontscan=treescan_north, backscan=treescan_south)\n",
    "                analysis.analysis(octfile, name=sim_name+'_East&West', frontscan=treescan_east, backscan=treescan_west)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "125ded6c",
   "metadata": {},
   "source": [
    "<a id='step5c'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc6a9aa6",
   "metadata": {},
   "source": [
    "#### Single simulation until MakeOct for Getting a PRETTY IMAGE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48d7c619",
   "metadata": {},
   "outputs": [],
   "source": [
    "tree_albedo = 0.165 # Wikipedia [0.15-0.18]\n",
    "trunk_x = 0.8 * ft2m\n",
    "trunk_y = trunk_x\n",
    "trunk_z = 1 * ft2m\n",
    "\n",
    "tree_x = 3 * ft2m\n",
    "tree_y = tree_x\n",
    "tree_z = 4 * ft2m\n",
    "\n",
    "\n",
    "clearance_height = clearance_heights[0]\n",
    "xgap = xgaps[-1]\n",
    "tilt = tilts[0]\n",
    "pitch = y * np.cos(np.radians(tilt))+Ds[-1]\n",
    "\n",
    "sim_name = (sim_general_name+'_ch_'+str(round(clearance_height,1))+\n",
    "                '_xgap_'+str(round(xgap,1))+\\\n",
    "                '_tilt_'+str(round(tilt,1))+\n",
    "                '_pitch_'+str(round(pitch,1)))\n",
    "\n",
    "\n",
    "demo = bifacial_radiance.RadianceObj(sim_name,str(testfolder))  \n",
    "demo.setGround(albedo)\n",
    "demo.readWeatherFile(epwfile)\n",
    "\n",
    "coffeeplant_x = (x+xgap)/2\n",
    "coffeeplant_y = pitch\n",
    "\n",
    "demo.gendaylit(4020)\n",
    "demo.makeModule(name=moduletype, x=x, y=y, xgap = xgap)\n",
    "sceneDict = {'tilt':tilt,'pitch':pitch,'clearance_height':clearance_height,'azimuth':azimuth, 'nMods': nMods, 'nRows': nRows} \n",
    "scene = demo.makeScene(moduletype=moduletype,sceneDict=sceneDict, hpc=hpc, radname = sim_name)\n",
    "\n",
    "\n",
    "for ii in range(0,3):\n",
    "    coffeeplant_x = (x+xgap)/2 + (x+xgap)*ii\n",
    "    for jj in range(0,3):\n",
    "        coffeeplant_y = pitch/2 + pitch*jj\n",
    "        name = 'tree'+str(ii)+str(jj)\n",
    "        text = '! genrev litesoil tube{}tree t*{} {} 32 | xform -t {} {} {}'.format('head'+str(ii)+str(jj),tree_z, tree_x/2.0, \n",
    "                                                                                         -trunk_x/2.0 + coffeeplant_x, \n",
    "                                                                                           -trunk_x/2.0 + coffeeplant_y, trunk_z)\n",
    "        text += '\\r\\n! genrev litesoil tube{}tree t*{} {} 32 | xform -t {} {} 0'.format('trunk'+str(ii)+str(jj),trunk_z, trunk_x/2.0, \n",
    "                                                                                         -trunk_x/2.0 + coffeeplant_x, \n",
    "                                                                                          -trunk_x/2.0 + coffeeplant_y)\n",
    "\n",
    "        customObject = demo.makeCustomObject(name,text)\n",
    "        demo.appendtoScene(radfile=scene.radfiles, customObject=customObject, text=\"!xform -rz 0\")\n",
    "\n",
    "octfile = demo.makeOct(octname = demo.basename , hpc=hpc)  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9501f465",
   "metadata": {},
   "source": [
    "#### Now you can view the Geometry by navigating on the terminal to the testfolder, and using the octfile name generated above\n",
    "\n",
    ">rvu -vf views\\front.vp -e .0265652 -vp 2 -21 2.5 -vd 0 1 0 Coffee_ch_1.8_xgap_1.2_tilt_18_pitch_2.2.oct"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daef615a",
   "metadata": {},
   "source": [
    "<a id='step6'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d42ba4a",
   "metadata": {},
   "source": [
    "## 6. Compile Results Trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f841c2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# irr_Coffee_ch_1.8_xgap_0.6_tilt_18_pitch_1.6_Front&Back.csv\n",
    "\n",
    "ch_all = []\n",
    "xgap_all = []\n",
    "tilt_all = []\n",
    "pitch_all = []\n",
    "NorthIrrad = []\n",
    "SouthIrrad = []\n",
    "EastIrrad = []\n",
    "WestIrrad = []\n",
    "\n",
    "\n",
    "ft2m = 0.3048\n",
    "clearance_heights = np.array([6.0, 8.0, 10.0])* ft2m\n",
    "xgaps = np.array([2, 3, 4]) * ft2m\n",
    "Ds = np.array([2, 3, 4]) * ft2m    # D is a variable that represents the spacing between rows, not-considering the collector areas.\n",
    "tilts = [18, 10]\n",
    "y = 1\n",
    "\n",
    "\n",
    "\n",
    "for ch in range (0, len(clearance_heights)):\n",
    "    \n",
    "    clearance_height = clearance_heights[ch]\n",
    "    for xx in range (0, len(xgaps)):\n",
    "        \n",
    "        xgap = xgaps[xx]\n",
    "\n",
    "        for tt in range (0, len(tilts)):\n",
    "        \n",
    "            tilt = tilts[tt]\n",
    "            for dd in range (0, len(Ds)):\n",
    "                pitch = y * np.cos(np.radians(tilt))+Ds[dd]\n",
    "\n",
    "                sim_name = ('irr_Coffee'+'_ch_'+str(round(clearance_height,1))+\n",
    "                                '_xgap_'+str(round(xgap,1))+\\\n",
    "                                '_tilt_'+str(round(tilt,1))+\n",
    "                                '_pitch_'+str(round(pitch,1))+'_North&South.csv')\n",
    "\n",
    "                sim_name2 = ('irr_Coffee'+'_ch_'+str(round(clearance_height,1))+\n",
    "                                '_xgap_'+str(round(xgap,1))+\\\n",
    "                                '_tilt_'+str(round(tilt,1))+\n",
    "                                '_pitch_'+str(round(pitch,1))+'_East&West.csv')\n",
    "\n",
    "                ch_all.append(clearance_height)\n",
    "                xgap_all.append(xgap)\n",
    "                tilt_all.append(tilt)\n",
    "                pitch_all.append(pitch)\n",
    "                data = pd.read_csv(os.path.join(resultsfolder, sim_name))\n",
    "                NorthIrrad.append(data['Wm2Front'].item())\n",
    "                SouthIrrad.append(data['Wm2Back'].item())\n",
    "                data = pd.read_csv(os.path.join(resultsfolder, sim_name2))\n",
    "                EastIrrad.append(data['Wm2Front'].item())\n",
    "                WestIrrad.append(data['Wm2Back'].item())\n",
    "\n",
    "\n",
    "ch_all = pd.Series(ch_all, name='clearance_height')\n",
    "xgap_all = pd.Series(xgap_all, name='xgap')\n",
    "tilt_all = pd.Series(tilt_all, name='tilt')\n",
    "pitch_all = pd.Series(pitch_all, name='pitch')\n",
    "NorthIrrad = pd.Series(NorthIrrad, name='NorthIrrad')\n",
    "SouthIrrad = pd.Series(SouthIrrad, name='SouthIrrad')\n",
    "EastIrrad = pd.Series(EastIrrad, name='EastIrrad')\n",
    "WestIrrad = pd.Series(WestIrrad, name='WestIrrad')\n",
    "\n",
    "df = pd.concat([ch_all, xgap_all, tilt_all, pitch_all, NorthIrrad, SouthIrrad, EastIrrad, WestIrrad], axis=1)\n",
    "df.to_csv(os.path.join(resultsfolder,'TREES.csv'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0423168b",
   "metadata": {},
   "outputs": [],
   "source": [
    "trees = pd.read_csv(os.path.join(resultsfolder, 'TREES.csv'))\n",
    "trees.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be73d173",
   "metadata": {},
   "outputs": [],
   "source": [
    "trees['TreeIrrad_percent_GHI'] = trees[['NorthIrrad','SouthIrrad','EastIrrad','WestIrrad']].mean(axis=1)*100/puerto_Rico_YEAR\n",
    "\n",
    "print(trees['TreeIrrad_percent_GHI'].min())\n",
    "print(trees['TreeIrrad_percent_GHI'].max())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8930dfc",
   "metadata": {},
   "source": [
    "<a id='step7'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44177c38",
   "metadata": {},
   "source": [
    "## 7. Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f2488a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "tilts_l = list(trees['tilt'].unique())\n",
    "ch_l = list(trees['clearance_height'].unique())\n",
    "print(tilts_l)\n",
    "print(ch_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9822688",
   "metadata": {},
   "outputs": [],
   "source": [
    "for tilt in tilts_l:\n",
    "    for clearance_height in ch_l:\n",
    "        df2=trees.loc[df['tilt']==tilts[1]]\n",
    "        df3 = df2.loc[df2['clearance_height']==clearance_heights[2]]\n",
    "        df3['pitch']=df3['pitch'].round(1)\n",
    "        df3['xgap']=df3['xgap'].round(1)\n",
    "\n",
    "        sns.set(font_scale=2) \n",
    "        table = df3.pivot('pitch', 'xgap', 'TreeIrrad_percent_GHI')\n",
    "        ax = sns.heatmap(table, cmap='hot', vmin = 22, vmax= 35, annot=True)\n",
    "        ax.invert_yaxis()\n",
    "        figtitle = 'Clearance Height ' + str(clearance_height/ft2m)+'  ft, Tilt ' + str(tilt) + '$^\\circ$'\n",
    "        plt.title(figtitle)\n",
    "        print(table)\n",
    "        plt.show()\n",
    "        \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}