8. More examples
[1]:
import openalea.plantgl.all as pgl_all
from lightvegemanager.LVM import LightVegeManager
from pgljupyter import SceneWidget
8.1. Sun positions
The print_sun function allows you to compare the sun position differences between RATP and CARIBU internal computation.
[38]:
from lightvegemanager.sun import print_sun
print(print_sun.__doc__)
Prints sun position ouputs from RATP and CARIBU algorithm with the same inputs
:param day: input day
:type day: int
:param hour: input hour
:type hour: int
:param coordinates: [latitude, longitude, timezone]
:type coordinates: list
:param truesolartime: activates true solar time or local time to compute sun position
:type truesolartime: bool
[11]:
print_sun(201, 12, [0.0, 0.0, 0.0], True)
print_sun(201, 12, [40.0, 0.0, 0.0], True)
print_sun(201, 16, [40.0, 12.0, 2.0], True)
print_sun(201, 16, [40.0, 12.0, 2.0], False)
print_sun(347, 16, [46.0, 0.0, 0.0], True)
--- SUN COORDONATES ---
--- Convention x+ = North, vector from sky to floor
--- azimut: south clockwise E = -90° W = 90° zenith: zenith = 0° horizon = 90°
--- day: 201 hour: 12 latitude: 0.00 °
--- true solar time
--- RATP ---
azimut: 180.000 zenith: 69.227
x: -0.354671 y: -0.000000 z: -0.934991
--- CARIBU ---
azimut: 180.000 zenith: 68.944
x: -0.359285 y: -0.000000 z: -0.933228
--- SUN COORDONATES ---
--- Convention x+ = North, vector from sky to floor
--- azimut: south clockwise E = -90° W = 90° zenith: zenith = 0° horizon = 90°
--- day: 201 hour: 12 latitude: 40.00 °
--- true solar time
--- RATP ---
azimut: 0.000 zenith: 70.773
x: 0.329307 y: -0.000000 z: -0.944223
--- CARIBU ---
azimut: 0.000 zenith: 71.148
x: 0.323132 y: -0.000000 z: -0.946354
--- SUN COORDONATES ---
--- Convention x+ = North, vector from sky to floor
--- azimut: south clockwise E = -90° W = 90° zenith: zenith = 0° horizon = 90°
--- day: 201 hour: 16 latitude: 40.00 °
--- true solar time
--- RATP ---
azimut: 87.963 zenith: 35.881
x: 0.028807 y: -0.809726 z: -0.586100
--- CARIBU ---
azimut: 88.295 zenith: 36.115
x: 0.024030 y: -0.807482 z: -0.589402
--- SUN COORDONATES ---
--- Convention x+ = North, vector from sky to floor
--- azimut: south clockwise E = -90° W = 90° zenith: zenith = 0° horizon = 90°
--- day: 201 hour: 16 latitude: 40.00 °
--- local time, true solar time is 18.70
--- RATP ---
azimut: 112.480 zenith: 5.641
x: -0.380514 y: -0.919537 z: -0.098291
--- CARIBU ---
azimut: 112.736 zenith: 5.895
x: -0.384436 y: -0.917420 z: -0.102706
--- SUN COORDONATES ---
--- Convention x+ = North, vector from sky to floor
--- azimut: south clockwise E = -90° W = 90° zenith: zenith = 0° horizon = 90°
--- day: 347 hour: 16 latitude: 46.00 °
--- true solar time
--- RATP ---
azimut: 52.853 zenith: 2.129
x: 0.603439 y: -0.796543 z: -0.037150
--- CARIBU ---
azimut: 53.169 zenith: 2.592
x: 0.598849 y: -0.799584 z: -0.045228
8.2. Two planes
A small example with two horizontal planes, in order to check if the xy orientation according to North-East-South-West.
[2]:
geom1 = pgl_all.FaceSet([(0,0,0),(2,0,0), (2,2,0),(0,2,0)],[range(4)]) # plaque dessous
geom2 = pgl_all.FaceSet([(0,1,1),(2,1,1), (2,3,1),(0,3,1)],[range(4)]) # plaque dessus, décalée en y (axe est-ouest) vers l'ouest
s = pgl_all.Scene([pgl_all.Shape(geom1, pgl_all.Material((250,0,0),1), 888),
pgl_all.Shape(geom2, pgl_all.Material((0,250,0),1), 999)])
[3]:
# input dict
geometry = {}
environment = {}
caribu_parameters = {}
# Paramètres pré-simulation
geometry["scenes"] = [s]
environment["coordinates"] = [0. ,0. ,0.] # latitude, longitude, timezone
environment["direct"] = True
environment["diffus"] = False
environment["reflected"] = False
environment["caribu opt"] = {}
environment["caribu opt"]["par"] = (0.10, ) # plaques opaques
environment["infinite"] = False
## Paramètres CARIBU ##
caribu_parameters["sun algo"] = "caribu"
caribu_parameters["caribu opt"] = {}
caribu_parameters["caribu opt"]["par"] = (0.10, ) # plaques opaques
# Déclaration de l'objet
lghtcaribu = LightVegeManager(environment=environment,
lightmodel="caribu",
lightmodel_parameters=caribu_parameters)
# création de la scène dans l'objet
lghtcaribu.build(geometry, global_scene_tesselate_level=5)
[39]:
# forçage arbitraire en µmol.m-2.s-1
PARi = 500
# jour arbitraire (21 septembre)
day = 264
# heure matin, soleil arrive de l'est, les 2 plaques reçoivent tout le rayonnement
hour = 16
print("--- day %i, hour %i"%(day, hour))
# Calcul du rayonnement
lghtcaribu.run(energy=PARi, day=day, hour=hour, parunit="micromol.m-2.s-1", truesolartime=True)
--- day 264, hour 16
[10]:
# visualisation
SceneWidget(lghtcaribu.to_plantGL(lighting=True),
position=(0.0, 0.0, 0.0),
size_display=(600, 400),
plane=False,
size_world = 10,
axes_helper=True)
[10]:
8.3. Examples with RATP and stem elements
Stems processing checking
-------------------------
Scene for verification is a set of horizontal planes in one voxel as so:
-------------
| ---- |
| ____ |
| ---- |
-------------
We compare the PAR mean on all the planes between CARIBU and RATP
Stems processing
* CARIBU
Computing of incident PAR is done with no transmittance, triangles are opaques
* RATP
triangles area are divided by 2
[34]:
geom1 = pgl_all.FaceSet([(0.1,0.1,0.96),(0.1,0.6,0.96), (0.6,0.6,0.96),(0.6,0.1,0.96)],[range(4)])
geom2 = pgl_all.FaceSet([(0.1,0.5,0.806),(0.1,0.9,0.806), (0.6,0.9,0.806),(0.6,0.5,0.806)],[range(4)])
geom3 = pgl_all.FaceSet([(0.5,0.1,0.646),(0.5,0.7,0.646), (0.9,0.7,0.646),(0.9,0.1,0.646)],[range(4)])
geom4 = pgl_all.FaceSet([(0.3,0.5,0.486),(0.3,0.9,0.486), (0.8,0.9,0.486),(0.8,0.5,0.486)],[range(4)])
geom5 = pgl_all.FaceSet([(0.3,0.1,0.32),(0.3,0.5,0.32), (0.8,0.5,0.32),(0.8,0.1,0.32)],[range(4)])
geom6 = pgl_all.FaceSet([(0.2,0.4,0.16),(0.2,0.9,0.16), (0.6,0.9,0.16),(0.6,0.4,0.16)],[range(4)])
s = pgl_all.Scene([pgl_all.Shape(geom1, pgl_all.Material((250,0,0),1), 888),
pgl_all.Shape(geom2, pgl_all.Material((250,0,0),1), 888),
pgl_all.Shape(geom3, pgl_all.Material((250,0,0),1), 888),
pgl_all.Shape(geom4, pgl_all.Material((250,0,0),1), 888),
pgl_all.Shape(geom5, pgl_all.Material((250,0,0),1), 888),
pgl_all.Shape(geom6, pgl_all.Material((250,0,0),1), 888)])
[35]:
geometry = {}
environment = {}
ratp_parameters = {}
caribu_parameters = {}
# Paramètres pré-simulation
geometry["scenes"] = [s]
geometry["stems id"] = [(888, 0)]
environment["coordinates"] = [0., 0., 0.] # latitude, longitude, timezone
environment["sky"] = "turtle46" # turtle à 46 directions par défaut
environment["diffus"] = False
environment["direct"] = True
environment["reflected"] = False
environment["infinite"] = False
## Paramètres CARIBU ##
caribu_parameters["sun algo"] = "caribu"
caribu_parameters["caribu opt"] = {}
caribu_parameters["caribu opt"]["par"] = (0.10, 0.05)
lghtcaribu = LightVegeManager(environment=environment,
lightmodel="caribu",
lightmodel_parameters=caribu_parameters)
lghtcaribu.build(geometry)
[36]:
## Paramètres RATP ##
dv = 1. # m
dx, dy, dz = dv, dv, dv # m
ratp_parameters["voxel size"] = [dx, dy, dz]
ratp_parameters["soil reflectance"] = [0., 0.]
ratp_parameters["mu"] = [1.]
ratp_parameters["tesselation level"] = 0
ratp_parameters["angle distrib algo"] = "compute global"
ratp_parameters["nb angle classes"] = 30
ratp_parameters["reflectance coefficients"] = [[0.1, 0.05]]
lghtratp = LightVegeManager(environment=environment,
lightmodel="ratp",
lightmodel_parameters=ratp_parameters)
lghtratp.build(geometry)
[28]:
# visualisation
SceneWidget(lghtratp.to_plantGL(printvoxels=True),
position=(0.0, 0.0, 0.0),
size_display=(600, 400),
plane=False,
size_world = 2,
axes_helper=True)
[28]:
[29]:
PARi=500
day=100.
hour=12
[ ]:
PARi=500
day=100.
hour=17.5
[31]:
lghtcaribu.run(energy=PARi, day=day, hour=hour, parunit="micromol.m-2.s-1", truesolartime=True)
print("=== CARIBU ===")
print(lghtcaribu.elements_outputs)
# visualisation
SceneWidget(lghtcaribu.to_plantGL(lighting=True),
position=(0.0, 0.0, 0.0),
size_display=(600, 400),
plane=False,
size_world = 2,
axes_helper=True)
=== CARIBU ===
Day Hour Organ VegetationType Area par Eabs par Ei
0 100.0 12 888 0 1.29 226.615845 251.795383
[31]:
[33]:
lghtratp.run(energy=PARi, day=day, hour=hour, parunit="micromol.m-2.s-1", truesolartime=True)
print("=== RATP ===")
print(lghtratp.elements_outputs)
# visualisation
SceneWidget(lghtratp.plantGL(lighting=True, printvoxels=True),
position=(0.0, 0.0, 0.0),
size_display=(600, 400),
plane=False,
size_world = 2,
axes_helper=True)
=== RATP ===
Day Hour Organ VegetationType Area PARa Intercepted \
0 100.0 12.0 888 2 1.29 377.546967 0.487036
Transmitted SunlitPAR SunlitArea ShadedPAR ShadedArea
0 0.487036 499.828644 0.487203 0.0 0.157797
[33]: