Source code for py_wake.site.wasp_grid_site

from py_wake import np
import xarray as xr
import pickle
import os
import glob
from _collections import defaultdict
import re
import copy
from py_wake.site.distance import TerrainFollowingDistance
from py_wake.site.xrsite import XRSite
from numpy import newaxis as na


[docs]class WaspGridSite(XRSite): """Site with non-uniform (different wind at different locations, e.g. complex non-flat terrain) weibull distributed wind speed. Data obtained from WAsP grid files"""
[docs] def __init__(self, ds, distance=TerrainFollowingDistance(), mode='valid'): """ Parameters ---------- ds : xarray dataset as returned by load_wasp_grd distance : Distance object, optional Alternatives are StraightDistance and TerrainFollowingDistance2 mode : {'valid', 'extrapolate'}, optional if valid, terrain elevation outside grid area is NAN if extrapolate, the terrain elevation at the grid border is returned outside the grid area """ self.use_WS_bins = True ds = ds.rename(A="Weibull_A", k="Weibull_k", f="Sector_frequency", spd='Speedup', orog_trn='Turning', elev='Elevation', sec='wd', z='h') ds = ds.assign_coords(wd=(ds.wd - 1) * (360 / len(ds.wd))) ds = ds.isel(x=np.where(~np.all(np.isnan(ds.Elevation), 1))[0], y=np.where(~np.all(np.isnan(ds.Elevation), 0))[0]) super().__init__(ds, distance=distance)
def _local_wind(self, localWind, ws_bins=None): lw = super()._local_wind(localWind, ws_bins) # ti is assumed to be the turbulence intensity given by CFD # (expected value of TI at 15m/s). The Normal Turbulence model # is used to calculate TI at different wind speed, # see footnote 4 at page 24 of IEC 61400-1 (2005) if 'time' in lw.coords: lw['TI_ilk'] = self.interp(self.ds.ti15ms, lw) * (.75 + 3.8 / lw.ws)[na, :, na] else: lw['TI_ilk'] = self.interp(self.ds.ti15ms, lw) * (.75 + 3.8 / lw.ws) return lw @classmethod def from_wasp_grd(cls, path, globstr='*.grd', speedup_using_pickle=True, distance=TerrainFollowingDistance(), mode='valid'): ds = load_wasp_grd(path, globstr, speedup_using_pickle) return WaspGridSite(ds, distance, mode)
def load_wasp_grd(path, globstr='*.grd', speedup_using_pickle=True): ''' Reader for WAsP .grd resource grid files. Parameters ---------- path: str path to file or directory containing goldwind excel files globstr: str string that is used to glob files if path is a directory. ''' var_name_dict = { 'Flow inclination': 'flow_inc', 'Mean speed': 'ws_mean', 'Meso roughness': 'meso_rgh', 'Obstacles speed': 'obst_spd', 'Orographic speed': 'orog_spd', 'Orographic turn': 'orog_trn', 'Power density': 'power_density', 'RIX': 'rix', 'Roughness changes': 'rgh_change', 'Roughness speed': 'rgh_spd', 'Sector frequency': 'f', 'Turbulence intensity': 'ti15ms', 'Weibull-A': 'A', 'Weibull-k': 'k', 'Elevation': 'elev', 'AEP': 'aep'} def _read_grd(filename): def _parse_line_floats(f): return [float(i) for i in f.readline().strip().split()] def _parse_line_ints(f): return [int(i) for i in f.readline().strip().split()] with open(filename, 'rb') as f: file_id = f.readline().strip().decode() nx, ny = _parse_line_ints(f) xl, xu = _parse_line_floats(f) yl, yu = _parse_line_floats(f) zl, zu = _parse_line_floats(f) values = np.array([l.split() for l in f.readlines() if l.strip() != b""], dtype=float) # around 8 times faster xarr = np.linspace(xl, xu, nx) yarr = np.linspace(yl, yu, ny) # note that the indexing of WAsP grd file is 'xy' type, i.e., # values.shape == (xarr.shape[0], yarr.shape[0]) # we need to transpose values to match the 'ij' indexing values = values.T ############# # note WAsP denotes unavailable values using very large numbers, here # we change them into np.nan, to avoid strange results. values[values > 1e20] = np.nan return xarr, yarr, values if speedup_using_pickle: if os.path.isdir(path): pkl_fn = os.path.join(path, os.path.split(os.path.dirname(path))[1] + '.pkl') if os.path.isfile(pkl_fn): try: with open(pkl_fn, 'rb') as f: return pickle.load(f) except Exception: print('loading %s failed. Loading from grid files instead' % pkl_fn) ds = load_wasp_grd(path, globstr, speedup_using_pickle=False) with open(pkl_fn, 'wb') as f: pickle.dump(ds, f, protocol=-1) return ds else: raise NotImplementedError if os.path.isdir(path): files = sorted(glob.glob(os.path.join(path, globstr))) else: raise Exception('Path was not a directory...') file_height_dict = defaultdict(list) pattern = re.compile(r'Sector (\w+|\d+) \s+ Height (\d+\.?\d*)m \s+ ([a-zA-Z0-9- ]+)') for f in files: sector, height, var_name = re.findall(pattern, f)[0] # print(sector, height, var_name) name = var_name_dict.get(var_name, var_name) file_height_dict[height].append((f, sector, name)) elev_avail = False first = True for height, files_subset in file_height_dict.items(): first_at_height = True for file, sector, var_name in files_subset: xarr, yarr, values = _read_grd(file) if sector == 'All': # Only 'All' sector has the elevation files. # So here we make sure that, when the elevation file # is read, it gets the right (x,y) coords/dims. if var_name == 'elev': elev_avail = True elev_vals = values elev_coords = {'x': xarr, 'y': yarr} elev_dims = ('x', 'y') continue else: var_name += '_tot' coords = {'x': xarr, 'y': yarr, 'z': np.array([float(height)])} dims = ('x', 'y', 'z') da = xr.DataArray(values[..., np.newaxis], coords=coords, dims=dims) else: coords = {'x': xarr, 'y': yarr, 'z': np.array([float(height)]), 'sec': np.array([int(sector)])} dims = ('x', 'y', 'z', 'sec') da = xr.DataArray(values[..., np.newaxis, np.newaxis], coords=coords, dims=dims) if first_at_height: ds_tmp = xr.Dataset({var_name: da}) first_at_height = False else: ds_tmp = xr.merge([ds_tmp, xr.Dataset({var_name: da})]) if first: ds = ds_tmp first = False else: ds = xr.concat([ds, ds_tmp], dim='z') if elev_avail: ds['elev'] = xr.DataArray(elev_vals, coords=elev_coords, dims=elev_dims) ############ # Calculate the compund speed-up factor based on orog_spd, rgh_spd # and obst_spd spd = 1 for dr in ds.data_vars: if dr in ['orog_spd', 'obst_spd', 'rgh_spd']: # spd *= np.where(ds.data_vars[dr].data > 1e20, 1, ds.data_vars[dr].data) spd *= np.where(np.isnan(ds.data_vars[dr].data), 1, ds.data_vars[dr].data) ds['spd'] = copy.deepcopy(ds['orog_spd']) ds['spd'].data = spd ############# # change the frequency from per sector to per deg # ds['f'].data = ds['f'].data * len(ds['f']['sec'].data) / 360.0 # ############# # # note WAsP denotes unavailable values using very large numbers, here # # we change them into np.nan, to avoid strange results. # for var in ds.data_vars: # ds[var].data = np.where(ds[var].data > 1e20, np.nan, ds[var].data) # make sure coords along z is asending order ds = ds.loc[{'z': sorted(ds.coords['z'].values)}] ###################################################################### if 'ti15ms' in ds and np.mean(ds['ti15ms']) > 1: ds['ti15ms'] *= 0.01 return ds def main(): if __name__ == '__main__': from py_wake.examples.data.ParqueFicticio import ParqueFicticio_path import matplotlib.pyplot as plt site = WaspGridSite.from_wasp_grd(ParqueFicticio_path, speedup_using_pickle=False) x, y = site.ds.x.values, site.ds.y.values Y, X = np.meshgrid(y, x) # plot elevation ax1, ax2 = plt.subplots(1, 2)[1] site.ds.Elevation.plot(ax=ax1, levels=100) i, j = 15, 12 ax1.plot(X[:, j], Y[:, j], 'r') ax1.axis('equal') Z = site.elevation_interpolator(X[:, j], Y[:, j], mode='extrapolate') ax2.plot(X[:, j], Z, 'r') # plot wind speed ax1, ax2 = plt.subplots(1, 2)[1] WS = (site.ds.ws * site.ds.Speedup).sel(ws=10, wd=0).interp(h=70) WS.plot(ax=ax1, levels=100) i, j = 15, 12 ax1.plot(X[:, j], Y[:, j], 'r') ax1.axis('equal') ax2.plot(X[:, j], WS[:, j], 'r') # plot shear z = np.arange(35, 200, 1) u_z = site.local_wind(x=[x[i]] * len(z), y=[y[i]] * len(z), h=z, wd=[0], ws=[10]).WS_ilk[:, 0, 0] plt.figure() plt.plot(u_z, z) plt.xlabel('Wind speed [m/s]') plt.ylabel('Height [m]') plt.show() main()