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 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.

In [1]: import xarray as xr

In [2]: import matplotlib.pyplot as plt

In [3]: import windkit as wk

In [4]: import pywasp as pw

Topography I/O

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

Supported topography file formats

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:

  • windkit.read_raster_map()

  • windkit.raster_map_to_file().

Similarly, all vector maps are read and written using:

  • windkit.read_vector_map()

  • 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:

In [5]: elev_map_raster = wk.read_raster_map(
   ...:     "source/tutorials/data/elev.tif",
   ...:     map_type="elevation"
   ...: )
   ...: 

In [6]: print(elev_map_raster)
<xarray.DataArray 'elevation' (south_north: 800, west_east: 800)>
[640000 values with dtype=float64]
Coordinates:
  * west_east    (west_east) float64 7.818 7.818 7.819 ... 8.482 8.483 8.483
  * south_north  (south_north) float64 56.11 56.11 56.11 ... 56.77 56.77 56.77
    crs          int64 0
Attributes:
    description:    height above sea level of a point on the surface of an or...
    long_name:      Terrain Elevation
    standard_name:  ground_level_altitude
    units:          m
    grid_mapping:   crs
    history:        2023-07-26T08:48:09:\twindkit==0.6.3\twk.read_raster_map(

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

In [7]: elev_map_raster.plot(cmap="viridis", vmin=0.0, vmax=60.0);
../_images/elev_raster_plot_example_ll.png

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

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

In [9]: elev_map_raster.plot(cmap="viridis", vmin=0.0, vmax=60.0);
../_images/elev_raster_plot_example_utm.png

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:

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

In [11]: elev_map_raster = wk.spatial.clip(elev_map_raster, bbox)

In [12]: elev_map_raster.plot(cmap="viridis", vmin=0.0, vmax=60.0);
../_images/elev_raster_plot_example_clipped.png

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 pywasp.io.rastermap_to_vectormap() function.

In [13]: elev_map_vector = pw.io.rastermap_to_vectormap(elev_map_raster)

In [14]: fig, ax = plt.subplots(1, 1, figsize=(6, 6))

In [15]: elev_map_raster.plot(ax=ax, cmap="Greys_r", vmin=0.0, vmax=60.0, add_colorbar=False)

In [16]: elev_map_vector.plot("elev", ax=ax, legend=True, legend_kwds={"label": "Elevation [m]"});
../_images/elev_vector_plot_example1.png

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.

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

In [18]: print(lc_map_vector)
                                              geometry  id_left  id_right
0    LINESTRING (520129.500 4608000.000, 520157.100...        0         0
1    LINESTRING (522602.400 4608000.000, 522613.500...        1         0
2    LINESTRING (518024.600 4608000.000, 518062.300...        0         1
3    LINESTRING (518062.300 4608074.500, 518160.500...        1         1
4    LINESTRING (522376.200 4608000.000, 522411.700...        2         1
..                                                 ...      ...       ...
889  LINESTRING (524303.300 4637000.000, 524300.400...        0         1
890  LINESTRING (519338.000 4636768.000, 519331.300...        3         0
891  LINESTRING (510688.600 4635512.000, 510675.500...        4         9
892  LINESTRING (510166.200 4637000.000, 510167.200...        1         9
893  LINESTRING (526903.900 4637000.000, 526904.100...        1         1

[894 rows x 3 columns]

Vector maps are geopandas.GeoDataFrame objects with 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 shapely.LineString’s’) and the elevation column “elev”.

Since vector maps are 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):

In [19]: fig, ax = plt.subplots(1, 1, figsize=(6, 6))

In [20]: ax.set_facecolor("black")

In [21]: cats = [3, 4]

In [22]: mask = lc_map_vector.id_left.isin(cats) | lc_map_vector.id_right.isin(cats)

In [23]: lc_map_vector.loc[mask].plot(
   ....:     "id_left",
   ....:     cmap="tab20",
   ....:     vmin=-0.5,
   ....:     vmax=19.5,
   ....:     linewidth=2,
   ....:     ax=ax
   ....: );
   ....: 

In [24]: 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))},
   ....: );
   ....:
../_images/elev_vector_plot_example2.png

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:

In [25]: print(lc_tbl)
id = LandCoverID
z0 = Roughness length (m)
d = Displacement height (m)
desc = Description
id	z0	d	desc
0	0.0300	0.0	
1	0.6000	0.0	
2	0.7000	0.0	
3	0.9000	0.0	
4	0.0000	0.0	
5	0.1000	0.0	
6	1.5000	0.0	
7	0.0100	0.0	
8	0.5000	0.0	
9	0.0500	0.0	
10	0.3000	0.0	
11	0.2000	0.0	
12	0.8000	0.0

Preparing Input for Terrain Modeling

To predict site effects from the Topography data with PyWAsP, the data is bundled together in the 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.

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

In [27]: topo_map = pw.wasp.TopographyMap(
   ....:     elev_map_vector,
   ....:     lc_map_vector,
   ....:     lc_tbl
   ....: )
   ....: 

In [28]: print(topo_map)
Roughness map

              PyWAsP map
              Contains 894 contour lines
              It has the extent of Bounds: (509764.0, 4608000.0) (529000.0, 4637000.0)
CRS: PROJCRS["WGS 84 / UTM zone 29N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 29N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-9,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["Between 12°W and 6°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Côte D'Ivoire (Ivory Coast). Faroe Islands. Guinea. Ireland. Jan Mayen. Mali. Mauritania. Morocco. Portugal. Sierra Leone. Spain. United Kingdom (UK). Western Sahara."],BBOX[0,-12,84,-6]],ID["EPSG",32629]]
              
Elevation map

              PyWAsP map
              Contains 582 contour lines
              It has the extent of Bounds: (509752.7, 4608000.0) (529000.0, 4637000.0)
CRS: PROJCRS["WGS 84 / UTM zone 29N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 29N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-9,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["Between 12°W and 6°W, northern hemisphere between equator and 84°N, onshore and offshore. Algeria. Côte D'Ivoire (Ivory Coast). Faroe Islands. Guinea. Ireland. Jan Mayen. Mali. Mauritania. Morocco. Portugal. Sierra Leone. Spain. United Kingdom (UK). Western Sahara."],BBOX[0,-12,84,-6]],ID["EPSG",32629]]

In 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.