{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Obstacle speedups\n\nThis example builds flat elevation and roughness maps from scratch, adds a\nbuilding with an interior courtyard and a porous row of trees, and plots the\nobstacle speedup over a resource grid.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Create simple maps\nFlat terrain with a single roughness length, so the only site effect in this\nexample is the shelter from the obstacles.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import geopandas as gpd\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport windkit as wk\nfrom shapely.geometry import Polygon, box\n\nimport pywasp as pw\n\n\ncrs = \"EPSG:32632\"\ncenter_x, center_y = 500_000.0, 6_200_000.0\nhalf_width = 1_000.0\nbbox = wk.spatial.BBox.from_cornerpts(\n    minx=center_x - half_width,\n    miny=center_y - half_width,\n    maxx=center_x + half_width,\n    maxy=center_y + half_width,\n    crs=crs,\n)\n\nelev_map = wk.create_vector_map(bbox, map_type=\"elevation\", elevation=0.0)\nroughness_map = wk.create_vector_map(\n    bbox, map_type=\"roughness\", roughness_change=(0.05, 0.05)\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Add a building and a porous tree row\nThe 60 m x 40 m building is solid and has an interior courtyard. The row of\ntrees south of it has the same 12 m height but 50% porosity, so the two\nwakes can be compared directly.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "building = Polygon(\n    box(center_x - 90, center_y + 20, center_x - 30, center_y + 60).exterior.coords,\n    holes=[\n        box(center_x - 70, center_y + 30, center_x - 50, center_y + 50).exterior.coords\n    ],\n)\ntree_row = box(center_x - 65, center_y - 110, center_x - 55, center_y - 10)\nobstacle_map = gpd.GeoDataFrame(\n    {\n        \"height\": [12.0, 12.0],\n        \"porosity\": [0.0, 0.5],\n    },\n    geometry=[building, tree_row],\n    crs=crs,\n)\n\ntopography = pw.wasp.TopographyMap(\n    elev_map,\n    roughness_map,\n    obstacle_map=obstacle_map,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Calculate and plot the obstacle speedup\nSite effects are directional; this plot shows the sector with wind from the\nwest at 6 m above ground, half the obstacle height. The solid building gives\na stronger and longer wake than the porous tree row, and the courtyard is\nsheltered by the surrounding building. Grid points inside the building and\nbelow its height have no obstacle speedup (NaN, with a warning), so they\nare blank. This ray- and sector-based view is consistent with the\nWAsP-shelter description in Pe\u00f1a et al. (2015), which discusses directional\nwake spreading and the reduction of shelter with increasing porosity. See\n[Shelter models and observations](https://backend.orbit.dtu.dk/ws/portalfiles/portal/122110497/Report_WP3_DTUversion.pdf).\n\nIncrease ``NR_POINTS_RESOURCE_GRID`` for a smoother picture; this also\nincreases the number of site-effect evaluations and therefore the run time.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "NR_POINTS_RESOURCE_GRID = 30\nHEIGHT = 6.0\ngrid = wk.spatial.create_cuboid(\n    west_east=np.linspace(center_x - 150, center_x + 250, NR_POINTS_RESOURCE_GRID),\n    south_north=np.linspace(center_y - 200, center_y + 200, NR_POINTS_RESOURCE_GRID),\n    height=[HEIGHT],\n    crs=crs,\n)\nsite_effects = topography.get_site_effects(grid, n_sectors=12)\nobstacle_speedup = site_effects[\"obstacle_speedups\"].sel(height=HEIGHT, sector=270)\n\nfig, ax = plt.subplots(figsize=(8, 6))\nobstacle_speedup.plot(ax=ax, x=\"west_east\", y=\"south_north\", cmap=\"RdYlBu_r\")\nobstacle_map.boundary.plot(ax=ax, color=\"black\", linewidth=1.5)\nax.set_title(f\"Obstacle speedup in sector 270\u00b0 at {HEIGHT:g} m a.g.l.\")\nax.set_aspect(\"equal\")\nplt.tight_layout()\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.13.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}