Statgen Calculation Service
The CalcStatgenService gives access to generalizing wind statistics from wind measurements. It is the first step towards a statistical based resource calculation, either using PARK or RESOURCE in windPRO. Likewise, it is the first step in windPRO scripting.
The service itself if relatively simple with the CalcStatgenService.GetStatgenCalc and CalcStatgenService.SetStatgenCalc for getting and setting calculations.
Additionally, it is possible to read wws (generalized wind statistic) files using StatgenService.GetWindstatDataFromFile.
Generating a generalized wind statistics in windPRO follows the following steps. The steps that are needed are the following:
Add or download terrain and roughness data
Make a statgen object for statistical wind resource calculations
Do a Statgen calculation to get the wind statistics
Add a statgen object for resource map calculation and link the wind statistics
Properties of the calculation object can be found here TApiCalcStatgen. The example below include generating a Statgen calculation from scratch. It is the same example as Resource Calculation Service and includes the generation of a resource map as well.
"""
Copyright 2023 EMD International
License for this script: MIT https://opensource.org/license/mit/
License for windPRO commercial software: https://www.emd-international.com/contact-us/general-terms-conditions-sale/
"""
import os
import numpy as np
from windproapi.utils import get_windpro_sample_path
from windproapi import WindProApi
from windproapi import nan_to_skipvalue
_windproapi = WindProApi()
_windproapi.start_windpro_random_port()
# This folder needs to exist on your machine
working_dir = os.path.join(get_windpro_sample_path('4.0'), 'Resources')
os.makedirs(working_dir, exist_ok=True)
project_path = os.path.join(working_dir, 'Resources.w36p')
# Get to all services
online_data_service = _windproapi.get_service('OnlineDataService')
project_service = _windproapi.get_service('ProjectService')
objects_service = _windproapi.get_service('ObjectsService')
calculation_service = _windproapi.get_service('CalculationService')
calc_statgen_service = _windproapi.get_service('CalcStatgenService')
calc_park_service = _windproapi.get_service('CalcParkService')
obj_site_data_service = _windproapi.get_service('ObjSiteDataService')
obj_meteo_service = _windproapi.get_service('ObjMeteoService')
obj_wtg_service = _windproapi.get_service('ObjWtgService')
calc_resource_service = _windproapi.get_service('CalcResourceService')
obj_wtg_areas_service = _windproapi.get_service('ObjWtgAreasService')
obj_elevation_grid_service = _windproapi.get_service('ObjElevationGridService')
windpro_service = _windproapi.get_service('WindproService')
factory = _windproapi.get_factory('WindproService')
# Position of the site center.
lng = 13
lat = 56.8
# Extend of the resgen around in m
east_west_resgen = 10_000
south_north_resgen = 10_000
# Heights used in the resource calculation
heights_resgen = [50, 100]
# resolution of the resource calculation
resolution = 100.0
# measurement height of the mast to be used for the resource calculation
measurement_height = 110.0
# Make a new project
project_service.NewProject(lng=lng, lat=lat, filename=project_path)
# Setting the project coordinate system to UTM
utm_epsg = project_service.GetUtmWgs84EpsgForProject()
project_service.SetProjectCoorEPSG(epsg=utm_epsg)
# Map data #
# Adding a layer for terrain data
terrain_layer = objects_service.AddLayer(layerName='Terrain data', parentFolderHandle=0)
objects_service.SetEditLayer(handle=terrain_layer)
# Getting elevation data sets
terrain_data_type = 'ODSTerrainGrid'
online_data_service.PrepareService(dataType=terrain_data_type, lat=lat, lon=lng)
terrain_services = online_data_service.GetServices(dataType=terrain_data_type)
print(terrain_services)
# Choosing a specific dataset and downloading the data
terrain_service_id = 'SWE50_Grid'
terrain_handle = online_data_service.DownloadHeightData(dataType=terrain_data_type,
implId=terrain_service_id,
userDesc='Terrain data grid',
lat=lat,
lon=lng,
width=30000,
height=30000,
useAsTin=True)
# Getting roughness dataset sets
roughness_data_type = 'ODSRoughnessGridToLine'
online_data_service.PrepareService(dataType=roughness_data_type, lat=lat, lon=lng)
roughness_service = online_data_service.GetServices(dataType=roughness_data_type)
print(roughness_service)
# Choosing a specific roughness dataset and downloading it
roughness_service_id = 'DataService_Rou_grid_GlobCover'
roughness_handle = online_data_service.DownloadRoughnessData(dataType=roughness_data_type,
implId=roughness_service_id,
userDesc='Roughness data',
lat=lat,
lon=lng,
width=30000,
height=30000)
# Adding forest layers
forest_layer = objects_service.AddLayer(layerName='Forest data', parentFolderHandle=0)
terrain_data_type = 'ODSObjHeightsGrid'
online_data_service.PrepareService(dataType=terrain_data_type, lat=lat, lon=lng)
available_data = online_data_service.GetServices(dataType=terrain_data_type)
handle_forest = online_data_service.DownloadHeightData(dataType=terrain_data_type,
implId='SLU_FOREST_TREEHEIGHTS_ServiceID',
userDesc='Forest Data',
lat=lat,
lon=lng,
width=30000,
height=30000,
useAsTin=False)
elevation_grid_obj = obj_elevation_grid_service.GetElevationGridObject(handle=handle_forest)
# Meteorological data #
# making a new layer for meteodata
meteo_layer = objects_service.AddLayer(layerName='Reference meteo data', parentFolderHandle=0)
objects_service.SetEditLayer(handle=terrain_layer)
# Getting climate data sets
# It is pretended that meso scale data can be used as measurement for simplicity of the script
meteo_data_type = 'ODSClimate'
online_data_service.PrepareService(dataType=meteo_data_type, lat=lat, lon=lng)
meteo_services = online_data_service.GetServices(dataType=meteo_data_type)
print(meteo_services)
# Selecting a specific dataset and downloading it
meteo_service_id = 'siEra5Basic'
meteo_handles = online_data_service.DownloadMeteoData(implId=meteo_service_id,
lat=lat,
lon=lng,
maxDist=30000,
numPoints=1,
fromYear=2014,
toYear=2015)
# Using the Meteo Service to get the
meteo_obj_era5 = obj_meteo_service.GetMeteoObject(handle=meteo_handles[0])
# Determining which height is closest to the selection height
meteo_height_diffs = []
meteo_heights = []
meteo_height_handles = []
for temp in meteo_obj_era5.Heights.TApiMeteoHeight:
meteo_height_diffs.append(temp.Height - measurement_height)
meteo_heights.append(temp.Height)
meteo_height_handles.append(temp.Handle)
index_min = np.argmin(abs(np.array(meteo_height_diffs)))
print('Closest height to {}m foud at {}m.'.format(measurement_height, meteo_heights[index_min]))
meteo_height_handle = meteo_height_handles[index_min]
# Flow modelling #
# Addign a layer of
site_data_layer = objects_service.AddLayer(layerName='Site data objects', parentFolderHandle=0)
objects_service.SetEditLayer(handle=site_data_layer)
# Adding a site data object
obj_site_data_statgen = objects_service.AddObject(apiObjType='SiteData',
lat=lat,
lng=lng,
userDesc='STATGEN for resource')
obj_site_data_service.SetPurpose(handle=obj_site_data_statgen.Handle, purpose='SdPurpStatgen')
# Connecting elevation and roughness data and
obj_site_data_service.TerrainLinkElevationAndRoughnessLine(handle=obj_site_data_statgen.Handle,
elevHandle=terrain_handle,
rouLineHandle=roughness_handle)
# Getting an update of the changes that have been made to the object back to python
obj_site_data_statgen = obj_site_data_service.GetSiteDataObject(obj_site_data_statgen.Handle)
# Making wind statistics
# Creating a new statgen calculation calculation and starting to edit it
wws_file = os.path.join(working_dir, 'wind_stat.wws')
statgen_handle = calculation_service.CreateEmpty(calcType='CalcStatgen')
calculation_service.OpenCalcForEdit(handle=statgen_handle)
# Getting the statgen calculation object and adding information on
statgen = calc_statgen_service.GetStatgenCalc()
statgen.Filename = wws_file
statgen.Name = 'Statgen for test'
statgen.SiteDataHandle = obj_site_data_statgen.Handle
statgen.MeteoHandle = meteo_obj_era5.Handle
statgen.MeteoHeiHandle = meteo_height_handle
nan_to_skipvalue(statgen)
calc_statgen_service.SetStatgenCalc(statgen)
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=statgen_handle)
# Reading generalized wind statistics for illustration
wind_stats = calc_statgen_service.GetWindstatDataFromFile(wws_file)
# Adding a site data object
obj_site_data_resource = objects_service.AddObject(apiObjType='SiteData',
lat=lat,
lng=lng,
userDesc='STATGEN for resource')
obj_site_data_service.SetPurpose(handle=obj_site_data_resource.Handle, purpose='SdPurpResource')
# Connecting elevation and roughness data and
obj_site_data_service.TerrainLinkElevationAndRoughnessLine(handle=obj_site_data_resource.Handle,
elevHandle=terrain_handle,
rouLineHandle=roughness_handle)
# Getting an update of the changes that have been made to the object back to python
obj_site_data_resource = obj_site_data_service.GetSiteDataObject(obj_site_data_resource.Handle)
# Setting wind statistigcs
wind_stats_files = factory.TApiSdWindstatFiles()
wind_stats_file = factory.TApiSdWindStatFile()
wind_stats_file.Filename = wws_file
wind_stats_files.TApiSdWindStatFile.append(wind_stats_file)
obj_site_data_resource.WindStatFiles = wind_stats_files
# Setting data back
nan_to_skipvalue(obj_site_data_resource)
obj_site_data_service.SetSiteDataObject(obj_site_data_resource)
obj_site_data_resource = obj_site_data_service.GetSiteDataObject(obj_site_data_resource.Handle)
print(obj_site_data_resource)
# Wind turbine area #
wtg_area_layer = objects_service.AddLayer(layerName='WTG Areas', parentFolderHandle=0)
# Making new WTG area object and getting its parameters
obj = objects_service.AddObject(apiObjType='WTGareas',
lat=lat,
lng=lng,
userDesc='Area for resgen')
wtg_areas_object = obj_wtg_areas_service.GetWtgAreasObject(obj.Handle)
# Paramter holding multiple WTG areas in a formate understandable for windPRO
wtg_areas_object.AreaList = factory.TApiObjWtgAreaList()
# Coordinate transformations
utm_zone = project_service.GetUtmWgs84EpsgForProject()
utm_coord_center = windpro_service.ConvertCoorEPSG(x=lng,
y=lat,
inEpsg=4326,
toEpsg=utm_zone)
# Calculate the coordinates of the corners of the wtg area
coord_wtg_areas = []
temp = [[east_west_resgen / 2, south_north_resgen / 2],
[-1 * east_west_resgen / 2, south_north_resgen / 2],
[-1 * east_west_resgen / 2, -1 * south_north_resgen / 2],
[east_west_resgen / 2, -1 * south_north_resgen / 2]]
for coord in temp:
coord_wtg_areas.append(windpro_service.ConvertCoorEPSG(x=utm_coord_center.X + coord[0],
y=utm_coord_center.Y + coord[1],
inEpsg=utm_zone,
toEpsg=4326))
# Adding data to the wind turbine area object
wtg_area = factory.TApiObjWtgArea()
wtg_area.Name = 'resgen'
# For the style of the line. See tns:TApiBrushStyle
wtg_area.Hatching = 'ApiBsClear'
# This will be the input to the wtg_are.Points.
wtg_area.Points = factory.TApiGPs()
for coord in coord_wtg_areas:
newGP = factory.TApiGP()
newGP.Lat = coord.Y
newGP.Lng = coord.X
wtg_area.Points.TApiGP.append(newGP)
# We can add e.g. a buffer zone
wtg_area.BufferZone = 500.0
# Appending the WTG area to the WTG areas. There can in principle be more than one.
nan_to_skipvalue(wtg_area)
wtg_areas_object.AreaList.TApiObjWtgArea.append(wtg_area)
# Giving data back to windPRO
nan_to_skipvalue(wtg_areas_object)
obj_wtg_areas_service.SetWtgAreasObject(wtg_areas_object)
# Adding the resgen calculation #
resoureceHandle = calculation_service.CreateEmpty(calcType='CalcResource')
calculation_service.OpenCalcForEdit(handle=resoureceHandle)
resource_calc = calc_resource_service.GetResourceCalc()
resource_calc.Name = 'Automate Resource Calculation'
resource_calc.ResourceType = 'ResourceTypeSiteData'
resource_calc.SiteDataHandle = obj_site_data_resource.Handle
resource_calc.Resolution = resolution
# setting heights
hub_heights = factory.TApiDoubleArray()
for h in heights_resgen:
hub_heights.double.append(h)
resource_calc.HubHeights = hub_heights
# Setting areas. Need to be a special datatype.
wtg_obj_handles = factory.TApiObjectHandles()
wtg_obj_handles.int.append(wtg_areas_object.Handle)
resource_calc.WTGAreaHandles = wtg_obj_handles
# Setting resource calculation so that it is ready
resource_calc.Filename = os.path.join(working_dir, 'resource.rsf')
nan_to_skipvalue(resource_calc)
calc_resource_service.SetResourceCalc(resource_calc)
resource_calc = calc_resource_service.GetResourceCalc()
print(resource_calc)
# Calculate resoureces
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=resource_calc.Handle)