Shadow Calculation Service
Warning
The display of result layers in shadow calculation are at present not working correctly. Manually recalculating will fix the issue.
CalcShadowService gives access to the shadow flickering module.
It only contains two function, CalcShadowService.GetShadowCalc for getting a shadow calculations and
CalcShadowService.SetShadowCalc for passing it back to windPRO.
With the general calculation service data and reports can be extracted from the shadow calculations.
Properties of the calculation object can be found here TApiCalcShadow.
Find below and example that gets data from a shadow calculation and exports the data (also including ZVI and noise calculations).
"""
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 import WindProApi, nan_to_skipvalue
from windproapi import compare_objects
from windproapi.utils import get_windpro_sample_path
from windproapi import compare_objects
def get_obj_from_list(objs, variable, value):
l_obj = []
for o in objs:
if o[variable] == value:
l_obj.append(o)
if len(l_obj) == 0:
raise ValueError(f'Could not find {value} in any {variable}.')
elif len(l_obj) == 1:
return l_obj[0]
else:
raise ValueError(f'More than one variable {variable} has value {value}.')
_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'), 'Environmental')
os.makedirs(working_dir, exist_ok=True)
project_path = os.path.join(working_dir, 'Environmental.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')
obj_wtg_service = _windproapi.get_service('ObjWtgService')
windpro_service = _windproapi.get_service('WindproService')
calc_ZVI_service = _windproapi.get_service('CalcZVIService')
wtg_explorer_service = _windproapi.get_service('WtgExplorerService')
obj_elevation_grid_service = _windproapi.get_service('ObjElevationGridService')
calc_decibel_service = _windproapi.get_service('CalcDecibelService')
obj_NSA_service = _windproapi.get_service('ObjNSAService')
calc_shadow_service = _windproapi.get_service('CalcShadowService')
factory = _windproapi.get_factory('WindproService')
# Position of the project
lng = 13
lat = 57
# Wind Farm locations in m from the site center
wtg_easting = [-100 + (i % 5) * 1000 for i in range(20)]
wtg_northing = [-4000 + np.floor(i / 5) * 1000 + (i % 5) * 100 for i in range(20)]
# 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
layer_terrain = objects_service.AddLayer(layerName='Terrain data', parentFolderHandle=0)
objects_service.SetEditLayer(handle=layer_terrain)
# 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=20000,
height=20000,
useAsTin=True)
# 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)
# Adding wind turbines
windfarm_layer = objects_service.AddLayer(layerName='Wind Farm', parentFolderHandle=0)
# Seaching for Vestas turbines. 0 indicate to take all not an actual filtering
wtgs = wtg_explorer_service.GetWtgsWithFilter(manufactor='VESTAS',
minRatedPower=0,
maxRatedPower=0,
minHubHeight=0,
maxHubHeight=0,
minRotorDiameter=0,
maxRotorDiameter=0)
wtg = [w for w in wtgs if 'VESTAS V90 2000 90.0 !O!' == w.DataName][0]
# Converting coordinate systems
utm_zone = project_service.GetUtmWgs84EpsgForProject()
# Converting from latitude/longitude using this for coordinate transformations
utm_center = windpro_service.ConvertCoorEPSG(x=lng,
y=lat,
inEpsg=4326,
toEpsg=utm_zone)
# Converted coordinates are equivalent to result from ConvertCoorEPSG
print('Easting: {:.1f}m, Northing: {:.1f}m'.format(utm_center.X, utm_center.Y))
for i, i_east, i_northing in zip(range(len(wtg_easting)), wtg_easting, wtg_northing):
position = windpro_service.ConvertCoorEPSG(x=utm_center.X + i_east,
y=utm_center.Y + i_northing,
inEpsg=utm_zone,
toEpsg=4326)
print('Longitude: {:.3f}, Latitude: {:.3f}'.format(position.X, position.Y))
new_obj = objects_service.AddObject(apiObjType='NewWTG',
lat=position.Y,
lng=position.X,
userDesc='New WTG obj')
wtgObj = obj_wtg_service.GetWtgObject(new_obj.Handle)
wtgObj.Filename = wtg.FileName
wtgObj.Hubheight = wtg.DefHubHeight
wtgObj.UserDescription = 'Turbine'
wtgObj.UserLabel = 'Turbine ' + str(i + 1).zfill(2)
nan_to_skipvalue(wtgObj)
res = obj_wtg_service.SetWtgObject(wtgObj)
wtgids = factory.TApiWtgIds()
for w in objects_service.GetObjects(apiObjType='NewWTG'):
dummy = factory.TApiWtgId()
dummy.Handle = w.Handle
dummy.Rowindex = 0
wtgids['TApiWtgId'].append(dummy)
# ZVI calculation
zvi_handle = calculation_service.CreateEmpty(calcType='CalcZVI')
# Open calculation and start setting properties
calculation_service.OpenCalcForEdit(handle=zvi_handle)
calc_zvi = calc_ZVI_service.GetZVICalc()
calc_zvi.Name = 'Script Generated'
# Changing the center of the calculation
calc_zvi.Center.Lat = lat + 0.01
calc_zvi.Center.Lng = lng
calc_zvi.CenterType = 'ApiZVICtUser'
# Setting the size of the calculation area
calc_zvi.Width = 20000
calc_zvi.Height = 20000
calc_zvi.Step = 250
# Adding all existing turbines to the calculation
# Adding the handle for the elevation grid and the forest height map that is to be used
grid_handles_int = factory.TApiObjectHandles()
grid_handles_int.int.append(terrain_handle)
grid_handles_int.int.append(handle_forest)
calc_zvi.GridHandles = grid_handles_int
# Adding all wind turbines
calc_zvi.NewWtgs = wtgids
nan_to_skipvalue(calc_zvi)
calc_ZVI_service.SetZVICalc(calc_zvi)
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=zvi_handle)
# Noise calculations
noise_layer = objects_service.AddLayer(layerName='Noise', parentFolderHandle=0)
# Adding a noise sensitive area that is the neighbouring city
lng_noise = 13.0637
lat_noise = 56.95564
new_obj = objects_service.AddObject(apiObjType='NSA',
lat=lat_noise,
lng=lng_noise,
userDesc='Script added')
nsa_obj = obj_NSA_service.GetNSAObject(new_obj.Handle)
# Adding more points to make it into an area
for lat, lng in zip([56.9637, 56.9582, 56.9417], [13.0680, 13.0891, 13.0733]):
point = factory.TApiGP()
point.Lat = lat
point.Lng = lng
nsa_obj.Points.TApiGP.append(point)
nsa_obj.FreeDefinable = True
nsa_obj.KindOfDemand = 'kodAbsolute'
nsa_obj.NoiseDemandNew = 43.0
nan_to_skipvalue(nsa_obj)
obj_NSA_service.SetNSAObject(nsa_obj)
# Compare this object to a newly generated object
# This is useful as there are many integers that define the behaviour of the NSA objects
new_obj = objects_service.AddObject(apiObjType='NSA',
lat=lat_noise+0.01,
lng=lng_noise,
userDesc='default script added')
nsa_obj_def = obj_NSA_service.GetNSAObject(new_obj.Handle)
nan_to_skipvalue(nsa_obj_def)
print(compare_objects(nsa_obj_def, nsa_obj))
objects_service.DeleteObject(new_obj.Handle)
# Noise calcualtion generated empty
noise_handle = calculation_service.CreateEmpty(calcType='CalcNoise')
# It is possible to change parameters in the calculation
calculation_service.OpenCalcForEdit(handle=noise_handle)
calc_noise = calc_decibel_service.GetDecibelCalc()
calc_noise.Name = 'Script generated'
# Adding all turbines turbines
calc_noise.NewWtgs = wtgids
# Adding all noise sensitice areas
nsas = factory.TApiWtgIds()
for obj in objects_service.GetObjects(apiObjType='NSA'):
nsa = factory.TApiWtgId()
nsa.Handle = obj.Handle
nsa.Rowindex = -1
nsas.TApiWtgId.append(nsa)
calc_noise.NSAs = nsas
nan_to_skipvalue(calc_noise)
calc_decibel_service.SetDecibelCalc(calc_noise)
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=noise_handle)
# When setting a calculation type many value of the calculation will automatically change
calculation_service.OpenCalcForEdit(handle=noise_handle)
calc_noise_edited = calc_decibel_service.GetDecibelCalc()
# It is possible to set another calculation type.
# This is only donw by numbers, so it is necessary to try it out in windPRO GUI and use that numeric
# 5 stands for 'Swedish, Jan 2002, Land'
calc_noise_edited.DecibelCalcType = 5
nan_to_skipvalue(calc_noise_edited)
calc_decibel_service.SetDecibelCalc(calc_noise_edited)
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=noise_handle)
# Shadow calculations
# Adding shadow receptors
shadow_layer = objects_service.AddLayer(layerName='Shadow', parentFolderHandle=0)
# Adding two shadow receptors representable for the two near by towns
# Keeping them to default settings fro simplicity
objects_service.AddObject(apiObjType='Shadow', lat=56.971983, lng=13.072742, userDesc='House 1')
objects_service.AddObject(apiObjType='Shadow', lat=56.965227, lng=13.063942, userDesc='House 2')
objects_service.AddObject(apiObjType='Shadow', lat=56.973484, lng=13.061792, userDesc='House 3')
# Creading shadow calcution
shadow_handle = calculation_service.CreateEmpty(calcType='CalcShadow')
calculation_service.OpenCalcForEdit(handle=shadow_handle)
calc_shadow = calc_shadow_service.GetShadowCalc()
# Add settings to a map as well
calc_shadow.Name = 'Script Generated'
calc_shadow.MapCalc = True
calc_shadow.WidthWest = 1000
calc_shadow.WidthEast = 7000
calc_shadow.HeightSouth = 7000
calc_shadow.HeightNorth = 1000
# Adding all shadow receptors to the calculation
receptors = factory.TApiWtgIds()
for obj in objects_service.GetObjects(apiObjType='Shadow'):
receptor = factory.TApiWtgId()
receptor.Handle = obj.Handle
receptor.Rowindex = -1
receptors.TApiWtgId.append(receptor)
calc_shadow.Receptors = receptors
# Adding all wtg objects
calc_shadow.NewWtgs = wtgids
nan_to_skipvalue(calc_shadow)
calc_shadow_service.SetShadowCalc(calc_shadow)
calculation_service.CloseCalc(save=True)
calculation_service.Calculate(handle=shadow_handle)