Commit 6ed1c9ce authored by Nikola Vasiljevic's avatar Nikola Vasiljevic
Browse files

fix issues with open vs read

parent 486fbb17
%% 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)
bwc = pw.binned_wind_climate.read_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:
 
> **__xarray note__** the `.where` method of `xarray.DataArray` leaves values fulfilling the condition as they are, while the values not fulfilling the condition is replaced by `np.nan`'s or the second argument of `.where`, if one is given (like above).
 
We will plot the roughness like we did the elevation. This time we do not need to tell `xarray` to use specific variables for `x` and `y`, instead we will let `xarray` use the defaults which happen to be `x` and `y` in the raster. Some additional code has been added to show the colorbar values in units of m instead of log(m).
 
%% Cell type:code id: tags:
 
``` python
# Plot log of surface roughness lengths
cs = rgh.plot(cmap='viridis',
vmin=-4.0,
vmax=0.2,
add_colorbar=False,
figsize=(14, 10))
 
# Add the colorbar ticks and values of the roughness length (z0) before log-transformation
cb = plt.colorbar(cs, ticks=[-4, -3, -2, -1, 0])
cb.ax.set_yticklabels([0.0001, 0.001, 0.01, 0.1, 1.0])
cb.set_label('z0')
 
# Add the Høvsøre mast location
plt.plot(loc_x, loc_y, 'or', ms=16)
 
# 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:
 
Again, we see the coastline clearly, and the water bodies south of the Høvsøre mast. To the east the plantation "Klosterhede plantage" is seen and north of it the city of "Lemvig".
 
## SITE EFFECTS
 
The terrain data are used to calculate site effects used to generalize the observed wind climate ("going up" in WAsP lingo) and to downscale to nearby locations ("going down").
 
The site effects are:
 
- Orographic speed-up (by sector)
- Internal boundar layer speed-up (by sector)
- Effective upstream surface roughness (by sector)
 
> **__pywasp note__** WAsP includes site effects due to obstacles, but we will not consider those here
 
### Create TopographyMap object
Let us investigate the site effects at the site. First we need to read in the raster data as `pywasp` `RasterMap` objects. This will allow us to convert the rasters to vectors, which is the most common data format used for terrain data in WAsP.
 
> **__pywasp note__** `pywasp` can work with raster maps directly, but it is recommended to convert them to vectors first
 
We will use `elev.tif` for the elevation instead of `elev.nc`, since the classmethod we will use `pywasp.RasterMap.from_file` currently expects a `.tif` file. The classmethod expects a `MapType` as the second argument, to know what kind of raster is being opened. Lastly, we will pass the proper epsg number as a keyword argument.
 
%% Cell type:code id: tags:
 
``` python
elev_map = pw.RasterMap.open_file('data/elev.tif', pw.MapType.Elevation, 4326)
rgh_map = pw.RasterMap.open_file('data/rough.tif', pw.MapType.Roughness, 3035)
elev_map = pw.RasterMap.read_file('data/elev.tif', pw.MapType.Elevation, 4326)
rgh_map = pw.RasterMap.read_file('data/rough.tif', pw.MapType.Roughness, 3035)
```
 
%% Cell type:markdown id: tags:
 
PyWAsP's `RasterMap` objects can be converted to `VectorMap` object using the `to_vectormap()` method. Notice that we "chain" the reprojection function to the elevation map vectorization to convert from the WGS84 Lat-Lon projection (4326) to the European Grid (3035).
 
> **pywasp note**: conversion to vector maps is a neccessary step for roughness data, and recommended when you need to reproject elevation data.
 
The elevation and roughness maps are the two components needed to make a `TopographyMap` object in `pywasp`, which is used to calculate the terrain effects.
 
%% Cell type:code id: tags:
 
``` python
elev_map = elev_map.to_vectormap().reproject(3035)
rgh_map = rgh_map.to_vectormap()
topo_map = pw.wasp.TopographyMap(elev_map, rgh_map)
```
 
%% Cell type:markdown id: tags:
 
### Roughness Rose
 
We will first explore the roughness elements upstream from the the Høvsøre mast for each of 12 sectors. We can do this by using `get_rou_rose_pt` method of `topo_map`, and pass in the location of the mast and the number of sectors we would like to use. It requires a WAsP configuration object (`pywasp.wasp.Config`) that stores the parameters. We will simply use the default values.
 
%% Cell type:code id: tags:
 
``` python
# calculate roughness elements from vetor using spider grid method, new WAsP method to be implemented
conf = pw.wasp.Config()
rou_rose = topo_map.get_rou_rose_pt(loc_x, loc_y, 12, conf,topo_map.get_elev_rose(loc_x, loc_y, 12, conf))
rou_rose.plot()
```
 
%% Output
 
 
%% Cell type:code id: tags:
 
``` python
# calculate roughness elements from vetor using slf method, classical WAsP method
conf = pw.wasp.Config()
conf.terrain.putpar(82, 0)
rou_rose = topo_map.get_rou_rose_pt(loc_x, loc_y, 12, conf)
rou_rose.plot()
```
 
%% Output