testscript.ipynb

{
    "cells": [
        {
            "attachments": {},
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Testscript.ipynb\n",
                "\n",
                "A Jupyter notebook for testing the `ArBuWe.py` module.\n",
                "Doubles as a *\"tutorial\"* of sorts."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Import module\n",
                "\n",
                "import ArBuWe as abw"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Initialize input parameters\n",
                "\n",
                "shapefile_path = \"data/FI.shp\"\n",
                "raster_path = \"data/gfa_res_curr_density_epsg4326.tif\"\n",
                "weather_start = \"2010-06-01\"\n",
                "weather_end = \"2010-06-02\""
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Try creating a `Shapefile` and inspect the uniform raster\n",
                "\n",
                "import matplotlib.pyplot as plt\n",
                "ec = (0,0,0,1) # Edge color for the vector plots (black)\n",
                "fc = (0,0,0,0) # Face color for the vector plots (transparent)\n",
                "lw = 0.2 # Line width for the vector plots\n",
                "\n",
                "shp = abw.Shapefile(shapefile_path)\n",
                "\n",
                "f = abw.plot_layout(shp, shp.raster, dpi=100, title=\"Shapefile\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Create weights for testing.\n",
                "\n",
                "#import random\n",
                "\n",
                "#weights = {'KU005':0.24, 'KU010':0.73, 'KU009':0.03,}\n",
                "#weights = {lid:random.uniform(0,1) for lid in shp.data.index}\n",
                "weights = {lid:1 for lid in shp.data.index}"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Try creating the cutout\n",
                "\n",
                "cutout = abw.prepare_cutout(shp, weather_start, weather_end)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Try preparing the `layout` raster and inspect nonzero data.\n",
                "\n",
                "raster, layout = abw.prepare_layout(shp, cutout, weights, raster_path, resampling=5)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Plot the raster and layout\n",
                "# This is commented out for the moment, as saving the layouts takes some time.\n",
                "\n",
                "#fig1 = abw.plot_layout(shp, raster, title=\"Raster\", dpi=300)\n",
                "#fig2 = abw.plot_layout(shp, layout, title=\"Layout\", dpi=300)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# ArBuWe development stuff\n",
                "\n",
                "The below cells were used for developing ArBuWe to be able to handle xarray-level initial heating/cooling demand calculations, as well as checking that the `atlite`\n",
                "tracking options work for irradiation functions as well.\n",
                "\n",
                "The below stuff is most likely useless, but I'm leaving it in as potential examples."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Raw diffuse irradiation\n",
                "\n",
                "diffi = cutout.irradiation(\n",
                "    orientation={\"slope\": 0.0, \"azimuth\": 0.0},\n",
                "    irradiation=\"diffuse\",\n",
                "    tracking=\"vertical\",\n",
                "    layout=layout, # Remove this line to get xarray level results.\n",
                ")\n",
                "diffi.plot()\n",
                "\n",
                "diffi_vertical = cutout.irradiation(\n",
                "    orientation={\"slope\": 90.0, \"azimuth\": 0.0},\n",
                "    irradiation=\"diffuse\",\n",
                "    tracking=\"vertical\",\n",
                "    layout=layout,\n",
                ")\n",
                "diffi_vertical.plot()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## DIFFUSE IRRADIATION IS LOWER FOR VERTICAL SURFACES!\n",
                "\n",
                "ArchetypeBuildingModel.jl doesn't account for this properly, causing more solar gains than there should be."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Raw ground irradiation?\n",
                "\n",
                "ground = cutout.irradiation(\n",
                "    orientation={\"slope\": 0.0, \"azimuth\": 0.0},\n",
                "    irradiation=\"ground\",\n",
                "    tracking=\"vertical\",\n",
                "    layout=layout, # Remove this line to get xarray level results.\n",
                ")\n",
                "ground.plot()\n",
                "\n",
                "ground_vertical = cutout.irradiation(\n",
                "    orientation={\"slope\": 90.0, \"azimuth\": 0.0}, # Ground irradiation is not affected by azimuth\n",
                "    irradiation=\"ground\",\n",
                "    tracking=\"vertical\",\n",
                "    layout=layout,\n",
                ")\n",
                "ground_vertical.plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Direct irradiation?\n",
                "\n",
                "dirs = {\n",
                "    \"horizontal\": (0.0, 0.0),\n",
                "    \"north\": (90.0, 0.0),\n",
                "    \"east\": (90.0, 90.0),\n",
                "    \"south\": (90.0, 180.0),\n",
                "    \"west\": (90.0, 270.0),\n",
                "}\n",
                "\n",
                "diri_none = {\n",
                "    dir: cutout.irradiation(\n",
                "        orientation={\"slope\": sl, \"azimuth\": az},\n",
                "        irradiation=\"direct\",\n",
                "        layout=layout, # Remove this line to get xarray level results.\n",
                "    )\n",
                "    for dir, (sl, az) in dirs.items()\n",
                "}\n",
                "\n",
                "diri_vertical = {\n",
                "    dir: cutout.irradiation(\n",
                "        orientation={\"slope\": sl, \"azimuth\": az},\n",
                "        irradiation=\"direct\",\n",
                "        tracking=\"vertical\",\n",
                "        layout=layout, # Remove this line to get xarray level results.\n",
                "    )\n",
                "    for dir, (sl, az) in dirs.items()\n",
                "}\n",
                "\n",
                "# Plot comparison, \"horizontal\" shouldn't depend on tracking.\n",
                "\n",
                "diri_none[\"horizontal\"].plot()\n",
                "diri_vertical[\"horizontal\"].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Plot comparison, no tracking should show 4 different cardinal directions (plus horizontal).\n",
                "\n",
                "for dir, (sl, az) in dirs.items():\n",
                "    diri_none[dir].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Plot comparison, vertical axis tracking should result in identical radiation for all cardinal directions (plus horizontal).\n",
                "\n",
                "for dir, (sl, az) in dirs.items():\n",
                "    diri_vertical[dir].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Plot comparison, vertical axis tracking vs no tracking?\n",
                "\n",
                "for dir, (sl, az) in dirs.items():\n",
                "    diri_none[dir].plot()\n",
                "    diri_vertical[dir].plot()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Seems reasonable enough.\n",
                "\n",
                "Vertical axis tracking seems to work as intended in `atlite`.\n",
                "Since we likely won't ever know the actual distribution of window and envelope surface\n",
                "areas towards the cardinal directions *(or any direction to be frank)* on the\n",
                "building-stock-scale, we might as well assume the vertical envelope areas to be more or less equally distributed.\n",
                "\n",
                "Thus, we should be able to use the `vertical` axis tracking option in `atlite` to\n",
                "calculate the effective solar irradiation on the vertical envelope.\n",
                "Naturally, this must take into account that not all of the envelope surface faces\n",
                "the correct way, but that shouldn't be an issue.\n",
                "\n",
                "\n",
                "### Next steps: Draft the calculations for the heating/cooling demand\n",
                "\n",
                "What coefficients do we actually need from ArBuMo, and how to efficiently manage\n",
                "the weather stuff for thermal mass nodes.\n",
                "\n",
                "\n",
                "#### Test preprocessing of weather\n",
                "\n",
                "First, we want to preprocess some weather data.\n",
                "Mainly, we're interested in the ambient air temperature as well as the effective total\n",
                "irradiation on horizontal and vertical surfaces."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Test weather data preprocessing.\n",
                "\n",
                "external_shading_coefficient = 0.7\n",
                "\n",
                "ambient_temperature, effective_irradiation = abw.preprocess_weather(\n",
                "    cutout, external_shading_coefficient\n",
                ")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "ambient_temperature[1].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "effective_irradiation['horizontal'][1].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "effective_irradiation['vertical'][1].plot()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "#### Generate test internal heat gains\n",
                "\n",
                "We need these for testing the initial heating demand calculations."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Define internal heat gains\n",
                "\n",
                "import numpy as np\n",
                "\n",
                "int_gains_W_m2 = np.array([2.17, 2.17, 2.17, 3.47, 5.25, 4.86, 3.20, 1.34, 0.79, 0.79, 0.79, 1.25, 3.29, 3.44, 3.90, 9.78, 8.60, 4.76, 12.75, 7.14, 4.21, 3.55, 2.17, 2.17])\n",
                "gross_floor_area_m2 = 135.56\n",
                "int_gains_W = int_gains_W_m2 * gross_floor_area_m2\n",
                "int_gains_W"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Expand the array to an xarray\n",
                "\n",
                "# Format to full-length timeseries\n",
                "time = ambient_temperature.time\n",
                "int_array = np.tile(int_gains_W, len(time)//len(int_gains_W))\n",
                "\n",
                "# Call the expansion function\n",
                "internal_heat_gains_W = abw.expand_to_xarray(\n",
                "    int_array,\n",
                "    ambient_temperature,\n",
                "    \"internal heat gains\",\n",
                "    \"W\"\n",
                ")\n",
                "internal_heat_gains_W"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "#### Test preprocessing of heating demand.\n",
                "\n",
                "So the next step is to actually calculate the initial heating demand on `xarray` level.\n",
                "Fortunately, with this type of an approach, the cooling demand is just the negative part of heating demand, so the same function can be used for both applications.\n",
                "\n",
                "NOTE! Internal heat gains will need some additional processing to fit the `xarray` format."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Define remaining required input data for the heating demand processing.\n",
                "\n",
                "self_discharge_coefficient_W_K = 0.0\n",
                "total_ambient_heat_transfer_coefficient_W_K = 118.453\n",
                "solar_heat_gain_convective_fraction = 0.6\n",
                "window_non_perpendicularity_correction_factor = 0.9\n",
                "total_normal_solar_energy_transmittance = 0.472\n",
                "vertical_window_surface_area_m2 = 25.295\n",
                "horizontal_window_surface_area_m2 = 0.0"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Define time-dependent set-point via xarray expansion.\n",
                "\n",
                "set_points = np.tile([284.15, 294.15], len(time)//2)\n",
                "set_point_K = abw.expand_to_xarray(\n",
                "    set_points,\n",
                "    ambient_temperature,\n",
                "    \"Heating set point\",\n",
                "    \"K\"\n",
                ")\n",
                "set_point_K"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Test processing initial heating demand.\n",
                "\n",
                "initial_demand = abw.process_initial_heating_demand(\n",
                "    set_point_K,\n",
                "    ambient_temperature,\n",
                "    effective_irradiation,\n",
                "    internal_heat_gains_W,\n",
                "    self_discharge_coefficient_W_K,\n",
                "    total_ambient_heat_transfer_coefficient_W_K,\n",
                "    solar_heat_gain_convective_fraction,\n",
                "    window_non_perpendicularity_correction_factor,\n",
                "    total_normal_solar_energy_transmittance,\n",
                "    vertical_window_surface_area_m2,\n",
                "    horizontal_window_surface_area_m2,\n",
                ")\n",
                "initial_demand[1].plot()"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "#### Test separation and aggregation of heating/cooling demand.\n",
                "\n",
                "Ok, seems like the heating/cooling demand calculations are working. Next, we'll have to see if we can actually separate aggregate the initial heating and cooling demands."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Separate initial demand into heating and cooling.\n",
                "\n",
                "import xarray as xr\n",
                "\n",
                "heating_demand = xr.where(initial_demand < 0.0, 0.0, initial_demand)\n",
                "cooling_demand = xr.where(initial_demand > 0.0, 0.0, -initial_demand)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Inspect initial demand\n",
                "\n",
                "initial_demand.plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Inspect heating demand\n",
                "\n",
                "heating_demand.plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cooling_demand.plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "initial_demand[2].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "heating_demand[2].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "cooling_demand[2].plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Try aggregating the weather?\n",
                "\n",
                "agg_ambient_temperature = abw.aggregate_xarray(ambient_temperature, layout)\n",
                "comp_temperature = cutout.temperature(layout=layout) + 273.15\n",
                "agg_ambient_temperature.plot(label=\"ArBuWe\")\n",
                "comp_temperature.plot(label=\"atlite\")\n",
                "plt.legend(loc=\"upper left\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Try aggregating the heating and cooling demands.\n",
                "\n",
                "agg_initial_demand = abw.aggregate_xarray(initial_demand, layout)\n",
                "agg_heating_demand = abw.aggregate_xarray(heating_demand, layout)\n",
                "agg_cooling_demand = abw.aggregate_xarray(cooling_demand, layout)"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "agg_initial_demand.plot(label=\"net\")\n",
                "agg_heating_demand.plot(label=\"heating\")\n",
                "agg_cooling_demand.plot(label=\"cooling\")\n",
                "plt.legend(loc=\"upper left\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Surprisingly, the aggregated heating/cooling demand calculations seem to be working as intended. The zig-zag is a consequence of time-varying heating/cooling set points, which one would never do in practise. This is purely for testing purposes.\n",
                "\n",
                "Furthermore, the aggregation of ambient temperature seems to result in identical values regardless of whether I use my own aggregation or PyPSA/atlite native aggregation functionality! *(Then again, it's basically just copy-paste from atlite, but at least I didn't mess anything up!)*\n",
                "\n",
                "#### Test the master processing function\n",
                "\n",
                "Finally, we'll need to test the master processing function aimed at the interface with ArBuMo.jl. Essentially, we'll need to return the initial heating and cooling demands, as well as the aggregated weather quantities for structural thermal mass dynamics in ArBuMo.jl."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Test the master processing function\n",
                "\n",
                "hd, cd, aT, irr = abw.aggregate_demand_and_weather(\n",
                "    shapefile_path,\n",
                "    weather_start,\n",
                "    weather_end,\n",
                "    weights,\n",
                "    external_shading_coefficient,\n",
                "    set_points,\n",
                "    set_points + 5,\n",
                "    int_array,\n",
                "    int_array - 5,\n",
                "    np.tile([self_discharge_coefficient_W_K], len(time)),\n",
                "    np.tile([total_ambient_heat_transfer_coefficient_W_K], len(time)),\n",
                "    np.tile([total_ambient_heat_transfer_coefficient_W_K * 2.0], len(time)),\n",
                "    solar_heat_gain_convective_fraction,\n",
                "    window_non_perpendicularity_correction_factor,\n",
                "    total_normal_solar_energy_transmittance,\n",
                "    vertical_window_surface_area_m2,\n",
                "    horizontal_window_surface_area_m2,\n",
                "    raster_path=None,\n",
                "    resampling=5,\n",
                "    filename=\"scope\",\n",
                "    save_layouts=True,\n",
                ")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "hd.plot(label=\"heating\")\n",
                "cd.plot(label=\"cooling\")\n",
                "plt.legend(loc=\"upper left\")"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "aT.plot()"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "irr[\"horizontal\"].plot(label=\"horizontal\")\n",
                "irr[\"vertical\"].plot(label=\"vertical\")\n",
                "plt.legend(loc=\"upper left\")"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "The horizontal irradiation appears much higher here compared to the vertical irradiation because the effective vertical irradiation is to be applied to all vertical surface area, not just the part that is actually facing the sun. For vertical surfaces matching the solar azimuth, the irradiation is Pi times higher.\n",
                "\n",
                "## Conclusion: it works?\n",
                "\n",
                "I almost can't believe it, but it seems that my initial heating and cooling demand processing is working as expected. I'm not sure if it will be fast enough to do what is required of it on the full Finnish/European scales, but you never know.\n",
                "There are actually things wrong with the old ArchetypeBuildingModel:\n",
                " - Diffuse irradiation is not identical for horizontal and vertical surfaces, which overestimates solar gains in ArchetypeBuildingModel.\n",
                " - Ground irradiation is neglegted in ArchetypeBuildingModel, which compensates for the above slightly, but not enough."
            ]
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3.8.14",
            "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.11.8"
        },
        "orig_nbformat": 4,
        "vscode": {
            "interpreter": {
                "hash": "64d5ff2237bad89141841d55ad00e98d09d1a3606a909d180b30c89cb53acad7"
            }
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}