Resource Calculation Service

The service CalcResourceService gives you access to generate resource maps in windPRO The service itself if relatively simple with the CalcResourceService.GetResourceCalc and CalcResourceService.SetResourceCalc for getting and setting calculations.

Generating a resource map from scratch involves the main part of the flow modelling in windPRO. 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

  • Make a WTG area to decide the limits of your calculation

  • Make a resource calculation linking the statgen for resource object and the WTG area object

Properties of the calculation object can be found here TApiCalcResource. It is more informative to understand this form an example. The code below illustrates how to make a resource map from scratch. it is the same code use to illustrate Statgen Calculation Service.

"""
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)