Commit 486fbb17 authored by Nikola Vasiljevic's avatar Nikola Vasiljevic
Browse files

fixed failing examples

parent e35cbc08
%% Cell type:markdown id: tags:
 
# <center> *PyWAsP* tutorial 1 <br><br> From observed wind climate to resource map <center>
 
<img align="right" src="data/updown.png">
 
A notebook by the PyWAsP team
 
## Introduction
 
This is an interactive tutorial that will introduce you to the basic WAsP workflow using PyWAsP.
 
After completing this tutorial you will be familiar with basic `pywasp` functionality, including:
 
- Opening and inspecting wind climate files.
- Opening and inspecting terrain data
- Calculating and inspecting site effects from the terrain data
- Creating a generalized wind climate from the observed wind climate
- Creating a resource map by downscaling the generalized wind climate
 
As you work your way through the notebook, make sure to run the python code in each cell in the order that they appear. You run the code by clicking on the cell (outlines around the cell should appear) and pressing `<shift> + <enter>` on the keyboard.
 
> **__notebook note__** if something looks wrong, or when errors occur, it can be helpful to restart the python kernel via the _kernel_ tab in the top
 
 
 
## Import packages
Usually the first step when writing a python program is importing standard and external packages we will need. For this analysis we will import `numpy`, `matplotlib`, `xarray`, `pywasp` itself, and some custom code made specificly for this notebook stored in a local file `pywasp_tutorial.py`.
 
- `numpy` is python's main numerical array package https://www.numpy.org/
- `matplotlib` is python's main plotting package https://matplotlib.org/
- `xarray` is a powerful high level package for labelled multi-dimensional arrays http://xarray.pydata.org
- `pywasp` docummentation is found at http://docs.wasp.dk/
- `pywasp_tutorila` contains a plotting function that we will use in this tutorial
 
%% Cell type:code id: tags:
 
``` python
import numpy as np
import matplotlib.pyplot as plt
 
import xarray as xr
 
import pywasp as pw
from pywasp_tutorial import plot_bwc
```
 
%% Cell type:markdown id: tags:
 
It is common to make a short alias of the package when importing using the `import - as -` syntax. The functions and classes of imported packages must be accessed through explicitly typing the package name or alias, e.g. `np.cos(3.14)` will use the cosine function from `numpy`.
 
Next, we will use a magic command that makes matplotlib figures appear in the notebook when a plot command is called.
 
%% Cell type:code id: tags:
 
``` python
%matplotlib inline
```
 
%% Cell type:markdown id: tags:
 
## Observed wind climate
 
Now we will import our (observed) binned wind climate from the Høvsøre mast. The data is stored in `hovsore.tab` in the `data` folder. It represents the observed wind climate for the period 2007-2015 at 100 m. above ground level.
 
> **__pywasp note__** `pywasp` includes functionality to both read and write wind climate files in many data formats, including ascii (`.tab`), xml (`.owc`), and netCDF (`.nc`).
 
The geospatial coordinate in the file is in the [European grid projection](https://en.wikipedia.org/wiki/European_grid#cite_note-1) ([EPSG:3035](https://epsg.io/3035)), which is the recommended projection for high resolution studies in Europe. We will tell `pywasp` this by explicitly passing a keyword argument `srs=3035` to the `open_bwc`.
 
%% Cell type:code id: tags:
 
``` python
bwc = pw.binned_wind_climate.open_bwc('data/hovsore.tab', srs=3035)
 
print(bwc)
```
 
%% Output
 
<xarray.Dataset>
Dimensions: (point3d: 1, sector: 12, wsbin: 30)
Coordinates:
* wsbin (wsbin) float64 0.5 1.5 2.5 3.5 4.5 ... 26.5 27.5 28.5 29.5
* sector (sector) float64 0.0 30.0 60.0 90.0 ... 240.0 270.0 300.0 330.0
south_north (point3d) float64 3.706e+06
west_east (point3d) float64 4.207e+06
height (point3d) float64 100.0
wsceil (wsbin) float64 1.0 2.0 3.0 4.0 5.0 ... 27.0 28.0 29.0 30.0
crs <U1 ''
Dimensions without coordinates: point3d
Data variables:
wdfreq (sector, point3d) float64 0.02869 0.04459 ... 0.09727 0.0156
wsfreq (wsbin, sector, point3d) float64 0.01167 0.01544 ... 0.0 0.0
amplif float64 1.0
offset float64 0.0
Attributes:
header: Høvsøre observed wind climate EPSG:3035
 
%% Cell type:markdown id: tags:
 
Notice that the `bwc` object is of type `<xarray.Dataset>` and that it contains four kinds of data:
 
1. **Dimensions**: core named dimensions
2. **Coordinates**: coordinate values along dimensions
3. **Data variables**: named arrays with data along 0-N named dimensions
4. **Attributes**: additional meta data attached to the dataset
 
> **__xarray note__** the primitive datatype and dimensions of each variable are also shown, along with a small sample of the data. `wsfreq` is a four-dimensional `float64` (double precision) variable along dimensions `(wsbin, sector, height, point)`
 
Xarray datasets wrap numpy arrays, annotating them with human-readable dimensions and coordinates, and allowing for easy subsetting, data manipulation, and plotting of the underlying data. An `xarray.Dataset` object is a collection of data variables, while each varible itself has type `xarray.DataArray`.
 
> **__xarray note__** Use the `.values` object attribute to access the underlying numpy array
 
Beyond the wind speed and wind direction distributions, the wind climate contains information about the height of the measurements (`height`) and the geospatial location (`west_east` and `south_north`), which in this case hold the location information in the projected coordinates of the EPSG:3035 projection.
 
In the cell below, we will store the location of the Høvsøre mast in variables `loc_x` and `loc_y` for later use.
 
%% Cell type:code id: tags:
 
``` python
loc_x = bwc['west_east']
loc_y = bwc['south_north']
```
 
%% Cell type:markdown id: tags:
 
The next step is to plot the wind rose and wind speed distributions in the binned wind climate. For convinience a plotting function `plot_bwc` has been implemented in `pywasp_tutorial` that will do this.
 
> **notebook note**: you can view the documentation for a function in jupyter notebooks by placing a `?` in front of the function, and you can get the entire function by using `??`.
 
%% Cell type:code id: tags:
 
``` python
plot_bwc(bwc);
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
As expected the prevailing wind direction at Høvsøre is from west, while the wind sector with the greatest relative wind speeds is the sector centered on `300.0` degrees and the lowest relative wind speeds comes from due north.
 
## Topography data
 
Now that we have inspected the wind climate, the next step is to load and inspect the terrain data.
 
### Terrain data (SRTM)
 
For this tutorial we have provided SRTM elevation data near Høvsøre in a netCDF file `elev.nc` in the `data` directory. We will use `xarray` to open this dataset with the `open_dataset` command. Since the dataset only contains one data variable `'elev'`, we will select that using the brackets syntax `['elev']`.
 
%% Cell type:code id: tags:
 
``` python
# Open rastermap of elevations
elev = xr.open_dataset('data/elev.nc')['elev']
 
print(elev)
```
 
%% Output
 
<xarray.DataArray 'elev' (band: 1, y: 800, x: 800)>
[640000 values with dtype=float64]
Coordinates:
* band (band) int64 1
* y (y) float64 56.11 56.11 56.11 56.11 ... 56.77 56.77 56.77 56.77
* x (x) float64 7.818 7.818 7.819 7.82 ... 8.481 8.482 8.483 8.483
x_proj (y, x) float64 ...
y_proj (y, x) float64 ...
Attributes:
transform: [ 8.33304385e-04 0.00000000e+00 7.81722527e+00 0.00000000...
crs: +init=epsg:4326
res: [0.0008333 0.0008333]
is_tiled: 0
nodatavals: -32768.0
 
%% Cell type:markdown id: tags:
 
Notice that while the `bwc` object was a collection of variables (`<xarray.Dataset>`), `elev` is a single data variable of type `xarray.DataArray`.
 
The elevation data is provided as a raster defined in the World Geodetic System 1984 (latitude, longitude) projection (EPSG:4326). For convinience, we have also provided the coordinates projected to the European grid projection (EPSG:3035) as the coordinates `x_proj` and `y_proj`.
 
To plot the elevation on a map, we will use the builtin `plot` method of `xarray.DataArray`. To tell `xarray` we want it to use the projected coordinates we pass in the keyword arguments `x='x_proj'` and `y='y_proj'`. Additionally, we tell it to use the `terrain` colormap, to limit the values between `-11` and `50` and to plot it on a canvas that is `14x10` inches large.
 
> **__xarray note__** `xarray.DataArray.plot` wraps matplotlib functions with sensible defaults, which makes it easy to inspect the data "on the go"
 
To highlight the location of the mast, we use the matplotlib `plot` function, and tell it to plot a red dot `ro` of size `16`.
 
To zoom in near the mast, extending 20 kilometers in each direction, we use matplotlib's `xlim` and `ylim` functions.
 
%% Cell type:code id: tags:
 
``` python
# Plot the elevation
 
elev.plot(x='x_proj',
y='y_proj',
cmap='terrain',
vmin=-11,
vmax=50,
figsize=(14, 10))
 
# Add the Høvsøre mast location to the plot
plt.plot(loc_x, loc_y, 'or', ms=16)
 
# Define the plot extent
plt_extent = 20000.0
 
# Set the plot extent
plt.xlim([loc_x - plt_extent, loc_x + plt_extent])
plt.ylim([loc_y - plt_extent, loc_y + plt_extent]);
```
 
%% Output
 
 
%% Cell type:markdown id: tags:
 
Denmark is flat - no surprise. But we can clearly see the coastline of Jutland stretching north-south.
 
### Roughness data (CORINE Landuse)
Next up is the surface roughness lengths. Here we start from the CORINE 2018 landuse dataset, which we then have translated to roughnesses using a lookup table (More about these lookup tables in Tutorial 3). The roughness map is provided as a GeoTiff file `rough.tif`, which is a common format in the GIS community. Conviniently, the CORINE dataset is using the European grid projection by default. Since we are using a GeoTiff file and not a NetCDF file, we need to open the raster using `open_rasterio` from `xarray`.
 
> **GeoTiff note**: GeoTiff files do not have multiple variables, therefore, `.tif` files open as `xarray.DataArray`' objects not `xarray.Dataset`'s, so we do not need to select the roughness data variable.
 
%% Cell type:code id: tags:
 
``` python
# Open rastermap of roughness lengths
rgh = xr.open_rasterio('data/rough.tif')
 
print(rgh)
```
 
%% Output
 
<xarray.DataArray (band: 1, y: 800, x: 800)>
[640000 values with dtype=float64]
Coordinates:
* band (band) int64 1
* y (y) float64 3.746e+06 3.746e+06 3.745e+06 ... 3.666e+06 3.666e+06
* x (x) float64 4.167e+06 4.167e+06 4.167e+06 ... 4.247e+06 4.247e+06
Attributes:
transform: (100.0, 0.0, 4166800.0, 0.0, -100.0, 3745700.0)
res: (100.0, 100.0)
is_tiled: 0
nodatavals: (nan,)
scales: (1.0,)
offsets: (0.0,)
 
%% Cell type:markdown id: tags:
 
Next we want to plot the roughnesses, but first, we will change the water roughness from 0 to 0.0001 and calculate the log, to make it easier to see differences in the map.
 
%% Cell type:code id: tags:
 
``` python
# Change water values from 0 to 0.0001 and take the log
rgh = rgh.where(rgh > 0.0, 0.0001)
rgh = np.log(rgh)
```
 
%% Cell type:markdown id: tags: