Engineering Forest Service

The engineering forest or displacement height model calculates effects from nearby forests. It can be used in Scalers, Statgen calculations or Park calculations. Via scripting you have access to engineering forest models and can add models EngForestService.AddEngForest(name) and modify existing models with the GetEngForests and setEngForest. An additional feature not directly available in windPRO is to perform displacement height calculations directly via CalculateDisplacementHeights.

The example below shows how to setup an engineering forest model using an online forest height map. An example is shown how so get the displacement heights by sector. Afterwards the engineering forest model used in: the setup of a Scaler (not used afterwards), a statistical wind resource calculation using WAsP, and a Park calcualtion with four turbines.

import os
import time
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()

# Makiing example
working_dir = os.path.join(get_windpro_sample_path('3.6'), 'ScalarService')
os.makedirs(working_dir, exist_ok=True)
project_path = os.path.join(working_dir, 'ScalarExample.w36p')

# Get to all services
scaler_service = _windproapi.get_service('ScalerService')
online_data_service = _windproapi.get_service('OnlineDataService')
project_service = _windproapi.get_service('ProjectService')
objects_service = _windproapi.get_service('ObjectsService')
obj_site_data_service = _windproapi.get_service('ObjSiteDataService')
eng_forest_service = _windproapi.get_service('EngForestService')
obj_elevation_grid_service = _windproapi.get_service('ObjElevationGridService')
obj_meteo_service = _windproapi.get_service('ObjMeteoService')
calc_park_service = _windproapi.get_service('CalcParkService')
calculation_service = _windproapi.get_service('CalculationService')
calc_statgen_service = _windproapi.get_service('CalcStatgenService')
obj_wtg_service = _windproapi.get_service('ObjWtgService')
factory = _windproapi.get_factory('WindproService')


def load_measurement_into_meteo(file, wls_f, m_obj):
    """
    Loading a single file with an already setup wls (inport filter) file into a meteobject.
    :param str file: path to the file containing data.
    :param str wls_f: Path to wls import filter file.
    :param zeep.objects.TApiWpObject m_obj: Meteo Object generated with 'ObjMeteoService'.
    :return:
    """
    # Meteo Import file needs to be set up
    meteo_import = factory.TApiMeteoImportSetup()
    meteo_import.LoggerSetupFile = wls_f
    meteo_import.DontUseInSetup = False
    meteo_import.IsOnlineData = False

    # For the import setup we need infomration on the folder structure and where the data is located.
    # Multiple folders could be included, in this case it is only one
    folder_items = factory.TApiFileFolderItems()
    folder_item = factory.TApiFileFolderItem()
    folder_item.FFType = factory.TFileFolderListType('fftFile')
    folder_item.FFName = file
    folder_item.FFMask = None
    folder_item.FFInclSubFolder = False
    folder_item.FFTreatFolderAsZip = False
    nan_to_skipvalue(folder_item)
    folder_items.TApiFileFolderItem.append(folder_item)

    # Add Folder structure to serach paths
    meteo_import.SearchPathes = folder_items

    # Add import setup to meteo object
    nan_to_skipvalue(meteo_import)
    obj_meteo_service.AddImportSetup(handle=m_obj.Handle, importSetup=meteo_import, loggerSetupFn=wls_f)

    # Automatically create heights. Needs to be done before loading data to determine what data should
    # be loaded in.
    obj_meteo_service.AutoCreateFromImportSetup(handle=m_obj.Handle)
    # Load the actual data
    obj_meteo_service.LoadAllData(handle=m_obj.Handle)


# Create a WindPRO project in a given location
lng = 13
lat = 57
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)

# 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=10000,
                                                       height=10000,
                                                       useAsTin=False)

elevation_grid_obj = obj_elevation_grid_service.GetElevationGridObject(handle=handle_forest)

# Adding the displacement height calculation
new_engforest_name = 'Displacement height from scripting'
new_eng_forest = eng_forest_service.AddEngForest(name=new_engforest_name)
new_eng_forest.HeightGridObjectHandle = handle_forest
new_eng_forest.InputType = 'fApihitHeightGrid'
new_eng_forest.HeightGridObjectLayerId = elevation_grid_obj.GridToHcSetup.Items.TApiGridToHCSetupItem[0].LayerID
nan_to_skipvalue(new_eng_forest)
eng_forest_service.SetEngForest(engForest=new_eng_forest)

# Displacement height calculations
N_sectors = 12
lats = [lat - 0.01, lat + 0.01, lat + 0.01, lat - 0.01]
lngs = [lng - 0.01, lng - 0.01, lng + 0.01, lng + 0.01]
calc_points = factory.TApiGPs()
for lat_i, lng_i in zip(lats, lngs):
    new_gp = factory.TApiGP()
    new_gp.Lat = lat_i
    new_gp.Lng = lng_i
    calc_points.TApiGP.append(new_gp)

# Interface to perform displacement height calculations directly.
# In the windPRO GUI this features needs a wtg at the location but is available in windPRO scripting
disp_heights = eng_forest_service.CalculateDisplacementHeights(engForestHandle=new_eng_forest.Handle,
                                                               siteDataHandle=0,
                                                               gps=calc_points,
                                                               sectors=N_sectors)
print('Displacement height calculations:')
for i, (lat, lng) in enumerate(zip(lats, lngs)):
    print(f'Location Latitude {lat}, Longitude {lng}:')
    print('Sector | displacement height [m]')
    print('--------------------------------')
    for j, dsp in enumerate(disp_heights[i].double):
        print('{:-6} |      {:2.2f}'.format(j, dsp))
    print('\n')

# Additional examples how to use the engineering forest model in Scalars, Statgen, or Park calculations
# Making a scalar that uses the engineering forest height
scaler = scaler_service.AddScaler(name='Engineering Forest Model')
scaler.DisplCalcHandle = new_eng_forest.Handle
scaler.ScalerType = 'ScalerTypeMastScalingGeo'
scaler.DisplType = 'ScalerDisplTypeCalc'
nan_to_skipvalue(scaler)
scaler_service.AddScaler(name=scaler)

# Using this scalar in a statical WAsP calculation with the displacement height.
# Adding mast data
mast_layer = objects_service.AddLayer(layerName='Mast data', parentFolderHandle=0)
objects_service.SetEditLayer(handle=mast_layer)
mast_file = os.path.join(os.path.dirname(__file__), 'data/meteo/New_Salem_2_South_60m_FINAL.csv')
wls_file = os.path.join(os.path.dirname(__file__), 'data/meteo/import_filter.wls')
new_obj = objects_service.AddObject(apiObjType='MeteoObjectData',
                                    lat=lat,
                                    lng=lng,
                                    userDesc='Mast from New Salem')
# Loading measurements with function defined above
load_measurement_into_meteo(mast_file, wls_file, new_obj)
meteo_object = obj_meteo_service.GetMeteoObject(handle=new_obj.Handle)

# Getting terrain data
terrain_layer = objects_service.AddLayer(layerName='Terrain data', parentFolderHandle=0)
objects_service.SetEditLayer(handle=terrain_layer)

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)
# 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=20000,
                                                        height=20000,
                                                        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)
# 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=20000,
                                                             height=20000)

# Site data and statgen calculations for 
site_data_layer = objects_service.AddLayer(layerName='Site Data Objects', parentFolderHandle=0)
objects_service.SetEditLayer(handle=site_data_layer)

# Statistical wind climate with a Stat Gen Calculation
obj_site_data_statgen = objects_service.AddObject(apiObjType='SiteData',
                                                  lat=lat,
                                                  lng=lng,
                                                  userDesc='STATGEN SiteData')
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)

# 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_object.Handle
statgen.MeteoHeiHandle = 1
# Adding displacement height calculation
statgen.DispHeightSetup.DispHeightType = 'DispHeightTypeDHC'
statgen.DispHeightSetup.DispHeightCalcHandle = new_eng_forest.Handle
nan_to_skipvalue(statgen)
calc_statgen_service.SetStatgenCalc(statgen)
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=statgen_handle)

# WAsP calculation 
siteDataWasp = objects_service.AddObject(apiObjType='SiteData',
                                         lat=lat,
                                         lng=lng,
                                         userDesc='WASP SiteData')
obj_site_data_service.SetPurpose(handle=siteDataWasp.Handle, purpose='SdPurpWasp')
obj_site_data_service.TerrainLinkElevationAndRoughnessLine(handle=siteDataWasp.Handle,
                                                           elevHandle=terrain_handle,
                                                           rouLineHandle=roughness_handle)
obj_site_data_service.SetWindstatFn(handle=siteDataWasp.Handle,
                                    fileName=wws_file)

# Adding wind turbines for PARK calculation
wtg_layer = objects_service.AddLayer(layerName='Wind Farm', parentFolderHandle=0)
objects_service.SetEditLayer(handle=wtg_layer)

wtgPath = os.path.join(os.path.dirname(__file__), 'data/SIEMENS SWT-2.3_test.wtg')

for i, i_lat, i_lng in zip(range(4), lats, lngs):
    new_obj = objects_service.AddObject(apiObjType='NewWTG',
                                        lat=i_lat,
                                        lng=i_lng,
                                        userDesc='New WTG obj')
    wtgObj = obj_wtg_service.GetWtgObject(new_obj.Handle)
    wtgObj.Filename = wtgPath
    wtgObj.Hubheight = 92.6
    wtgObj.UserDescription = 'Turbine'
    wtgObj.UserLabel = 'Turbine ' + str(i + 1)
    nan_to_skipvalue(wtgObj)
    res = obj_wtg_service.SetWtgObject(wtgObj)

# Adding a PARK calculation
# Creating a new PARK calculation and opening it for editing
park_handle = calculation_service.CreateEmpty(calcType='CalcPark')
calculation_service.OpenCalcForEdit(handle=park_handle)
calc_park = calc_park_service.GetParkCalc()
calc_park.ParkCalcType = 'ParkTypeSiteData'
calc_park.Name = 'Park for test'

# Adding displacement height calculations
calc_park.StatData.DispHeightSetup.DispHeightType = 'DispHeightTypeDHC'
calc_park.StatData.DispHeightSetup.DispHeightCalcHandle = new_eng_forest.Handle

# Adding WAsP handle
intArr = factory.TApiObjectHandles()
intArr['int'].append(siteDataWasp.Handle)
calc_park.StatData.SiteDataObjects = intArr

# Making an array of WTGs that can be added to the park
wtgids = factory.TApiWtgIds()
for w in objects_service.GetObjects(apiObjType='NewWTG'):
    dummy = calc_park_service.GetTApiWtgId()
    dummy.Handle = w.Handle
    dummy.Rowindex = 0
    wtgids['TApiWtgId'].append(dummy)

# Adding turbines
calc_park.NewWtgs = wtgids
calc_park.ExistWtgs = factory.TApiWtgIds()
nan_to_skipvalue(calc_park)

# Sending the PARK calculation back to windPRO and calculating
calc_park_service.SetParkCalc(calc_park)
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=park_handle)