.. _topography:

==========
Topography
==========

Topographic information is obviously an important input to WAsP's flow models.
Specifically, WAsP's terrain model takes as input maps of surface elevation,
surface roughness, and optionally displacement height.

In PyWAsP, instead of dealing with surface roughness maps and displacement height maps as completely separate layers,
we deal with them as land cover maps, where patches of land are associated with specific land cover classes,
each having specific surface characteristics (i.e. the surface roughness and displacement height).

Maps of surface roughness can be read by :doc:`WindKit <windkit:windkit>`, but are converted directly to land cover classes and an associated look-up table of
surface characteristics.

To start working with topographic data, we will import ``xarray``, ``matplotlib``, ``windkit``, and ``pywasp``.

.. ipython:: python

    import xarray as xr
    import matplotlib.pyplot as plt
    import windkit as wk
    import pywasp as pw


Topography I/O
========================

Reading and writing topography files is handled by :doc:`WindKit <windkit:windkit>`. Currently,
windkit supports the following file formats:

.. list-table:: Supported topography file formats
   :widths: 2, 3, 2, 2
   :align: center
   :header-rows: 1

   * - Format
     - Raster/Vector
     - Read
     - Write
   * - ``.nc``
     - Raster
     - yes
     - yes
   * - ``.grd``
     - Raster
     - yes
     - yes
   * - ``.tif``
     - Raster
     - yes
     - yes
   * - ``.map``
     - Vector
     - yes
     - yes
   * - ``.gml``
     - Vector
     - yes
     - yes
   * - ``.gpkg``
     - Vector
     - yes
     - yes

All raster maps, regardless of the format, are read and written using:
 - :py:func:`windkit.read_raster_map`
 - :py:func:`windkit.raster_map_to_file`.

Similarly, all vector maps are read and written using:
 - :py:func:`windkit.read_vector_map`
 - :py:func:`windkit.vector_map_to_file`


Raster Maps
========================

To start, let's read some raster elevation data. The data covers part of the west coast of Jutland, Denmark:

.. ipython:: python

    elev_map_raster = wk.read_raster_map(
        "source/tutorials/data/elev.tif",
        map_type="elevation"
    )
    print(elev_map_raster)

Since the data is a ``xr.DataArray`` it is easy to visualize:

.. ipython:: python

    @savefig elev_raster_plot_example_ll.png width=6in
    elev_map_raster.plot(cmap="viridis", vmin=0.0, vmax=60.0);

The raster is in latitude-longitude projection (EPSG:4326).
To reproject the data to a new raster in UTM coordinates, :py:func:`windkit.spatial.warp` can be used:

.. ipython:: python

    elev_map_raster = wk.spatial.warp(elev_map_raster, to_crs="EPSG:25832", method="cubic")

    @savefig elev_raster_plot_example_utm.png width=6in
    elev_map_raster.plot(cmap="viridis", vmin=0.0, vmax=60.0);

Warping to a new raster can create cells with missing data along the edges,
so to clip the raster to only include real data, we will create a bounding box and clip to its bounds:

.. ipython:: python

    bbox = wk.spatial.BBox.from_cornerpts(
        minx=433_000,
        miny=6_245_000.0,
        maxx=468_000,
        maxy=6_280_000,
        crs="EPSG:25832"
    )

    elev_map_raster = wk.spatial.clip(elev_map_raster, bbox)

    @savefig elev_raster_plot_example_clipped.png width=6in
    elev_map_raster.plot(cmap="viridis", vmin=0.0, vmax=60.0);


While WindKit and PyWAsP have great support for working with both raster and vector maps, PyWAsP works best with
vector maps as the input. Therefore, PyWAsP enables conversion from raster maps to vector maps, through the
:py:func:`pywasp.io.rastermap_to_vectormap` function.

.. ipython:: python

    elev_map_vector = pw.io.rastermap_to_vectormap(elev_map_raster)

    fig, ax = plt.subplots(1, 1, figsize=(6, 6))
    elev_map_raster.plot(ax=ax, cmap="Greys_r", vmin=0.0, vmax=60.0, add_colorbar=False)

    @savefig elev_vector_plot_example1.png width=6in
    elev_map_vector.plot("elev", ax=ax, legend=True, legend_kwds={"label": "Elevation [m]"});


Vector Maps
========================

To further illustrate how to work with vector maps, let's read some land cover data from the Serra Santa Luzia site.
Since the data comes in a ``.map``, which does not hold information about the
coordinate reference system, and, because it can hold both roughness and elevation information, we
have to provide the ``crs`` value and ask for the "roughness" part of the data via the ``map_type="roughness"`` argument.
Remember, the roughness map is converted to land cover classes by windkit.

.. ipython:: python

    lc_map_vector, lc_tbl = wk.read_vector_map(
        "source/tutorials/data/SerraSantaLuzia.map",
        map_type="roughness",
        crs="EPSG:32629",
    )
    print(lc_map_vector)

Vector maps are :py:class:`geopandas.GeoDataFrame` objects with :py:class:`shapely.LineString` geometries associated with either land cover boundaries or elevation levels.
Here, we have just two columns of data, the geometry column (for the :py:class:`shapely.LineString`'s') and the elevation column "elev".

Since vector maps are :py:class:`geopandas.GeoDataFrame` objects, they are easy to inspect and plot.
To avoid too much clutter, we will only plot lines where some specific classes are present (ocean and one forest category):

.. ipython:: python

    fig, ax = plt.subplots(1, 1, figsize=(6, 6))
    ax.set_facecolor("black")
    cats = [3, 4]
    mask = lc_map_vector.id_left.isin(cats) | lc_map_vector.id_right.isin(cats)
    lc_map_vector.loc[mask].plot(
        "id_left",
        cmap="tab20",
        vmin=-0.5,
        vmax=19.5,
        linewidth=2,
        ax=ax
    );
    @savefig elev_vector_plot_example2.png width=8in align=center
    lc_map_vector.loc[mask].plot(
        "id_right",
        cmap="tab20",
        vmin=-0.5,
        vmax=19.5,
        ax=ax,
        linewidth=1,
        linestyle="--",
        legend=True,
        legend_kwds={"label": "Land cover category", "ticks": list(range(13))},
    );

The way we plot the roughness-change lines, each line has two colors: the index for the "left" and "right" land cover category.
The lines represent jumps from one type of land cover to another.
The land cover table holds the information about the classes and their associated surface roughnesses and displacement heights:

.. ipython:: python

    print(lc_tbl)


Preparing Input for Terrain Modeling
====================================

To predict site effects from the Topography data with PyWAsP, the data is bundled together
in the :class:`pywasp.wasp.TopographyMap` PyWAsP class object.
It takes as input an elevation vector map, a land cover vector map, and a land cover table, as input.


.. ipython:: python

    elev_map_vector = wk.read_vector_map(
        "source/tutorials/data/SerraSantaLuzia.map",
        map_type="elevation",
        crs="EPSG:32629",
    )

    topo_map = pw.wasp.TopographyMap(
        elev_map_vector,
        lc_map_vector,
        lc_tbl
    )
    print(topo_map)


Topography maps can be saved to, and loaded from, a ZipFile archive using the :py:method:`pywasp.wasp.TopographyMap.save` and :py:func:`pywasp.wasp.TopographyMap.load` methods.

In :ref:`wasp_flow_model` we describe the WAsP models in more details and show how terrain effects can be calculated
from a ``topo_map``, like the one above.