diff --git a/src/alinea/pyratp/RatpScene.py b/src/alinea/pyratp/RatpScene.py index 7b2355d9..733f545d 100644 --- a/src/alinea/pyratp/RatpScene.py +++ b/src/alinea/pyratp/RatpScene.py @@ -431,7 +431,7 @@ def _dist(inc): return grid, vox_id, sh_id, s - def do_irradiation(self, rleaf=[0.1], rsoil=0.20, doy=1, hour=12, Rglob=1, Rdif=1, mu=None, sources=None): + def do_irradiation(self, rleaf=(0.1,), rsoil=0.20, doy=1, hour=12, Rglob=1, Rdif=1, mu=None, sources=None): """ Run a simulation of light interception for one wavelength Parameters: @@ -457,7 +457,7 @@ def do_irradiation(self, rleaf=[0.1], rsoil=0.20, doy=1, hour=12, Rglob=1, Rdif= vegetation = Vegetation.initialise(entities, nblomin=1) - if sources == None: + if sources is None: sky = Skyvault.initialise() else: el, az, w = sources @@ -474,31 +474,29 @@ def do_irradiation(self, rleaf=[0.1], rsoil=0.20, doy=1, hour=12, Rglob=1, Rdif= met = MicroMeteo.initialise(doy=doy, hour=hour, Rglob=Rglob, Rdif=Rdif) - res = runRATP.DoIrradiation(grid, vegetation, sky, met) + dfvox = runRATP.DoIrradiation(grid, vegetation, sky, met) - VegetationType,Iteration,day,hour,VoxelId,ShadedPAR,SunlitPAR,ShadedArea,SunlitArea= res.T - # 'PAR' is expected in Watt.m-2 in RATP input, whereas output is in micromol => convert back to W.m2 (cf shortwavebalance, line 306) - dfvox = pandas.DataFrame({'VegetationType':VegetationType, - 'Iteration':Iteration, - 'day':day, - 'hour':hour, - 'VoxelId':VoxelId, - 'ShadedPAR':ShadedPAR / 4.6, - 'SunlitPAR':SunlitPAR / 4.6, - 'ShadedArea':ShadedArea, - 'SunlitArea': SunlitArea, - 'Area': ShadedArea + SunlitArea, - 'PAR': (ShadedPAR * ShadedArea + SunlitPAR * SunlitArea) / (ShadedArea + SunlitArea) / 4.6, - }) - dfvox = dfvox[dfvox['VegetationType'] > 0] - dfvox = pandas.merge(dfvox, self.voxel_map.loc[:, - ('VoxelId', 'x', 'y', 'z', 'Volume')]) + # 'PAR' is expected in Watt.m-2 in RATP input, + # whereas output is in micromol + # => convert back to W.m2 (cf shortwavebalance, line 306) + dfvox['ShadedPAR'] /= 4.6 + dfvox['SunlitPAR'] /= 4.6 + + dfvox['Area'] = dfvox['ShadedArea'] + dfvox['SunlitArea'] + dfvox['PAR'] = (dfvox['ShadedPAR'] * dfvox['ShadedArea'] + + dfvox['SunlitPAR'] * dfvox['SunlitArea']) / dfvox['Area'] + + dfvox = dfvox.reset_index(drop=True) # for backward compatibility reasons + # TODO avoid and perform concat instead of merge + dfvox = pandas.merge(dfvox, self.voxel_map.loc[:, ('VoxelId', 'x', 'y', 'z', 'Volume')]) index = range(len(voxel_id)) dfmap = pandas.DataFrame( - {'primitive_index': index, 'shape_id': shape_id, + {'primitive_index': index, + 'shape_id': shape_id, 'VoxelId': voxel_id, 'VegetationType': [self.entity[sh_id] for sh_id in shape_id], - 'primitive_area': areas}) + 'primitive_area': areas} + ) return dfvox, dfmap diff --git a/src/alinea/pyratp/interface/ratp_scene.py b/src/alinea/pyratp/interface/ratp_scene.py index eee24d62..329d1859 100644 --- a/src/alinea/pyratp/interface/ratp_scene.py +++ b/src/alinea/pyratp/interface/ratp_scene.py @@ -368,28 +368,20 @@ def do_irradiation(self, sun_sources=None, sky_sources=None, azmoy = numpy.array(azmoy) - 180 + self.orientation sky = Skyvault.initialise(hmoy=hmoy, azmoy=azmoy, omega=omega, pc=pc) met = MicroMeteo.initialise(doy=1, hour=12, Rglob=ghi, Rdif=ghi) - res = runRATP.DoIrradiation(self.ratp_grid, vegetation, sky, met) - - VegetationType, Iteration, day, hour, VoxelId, ShadedPAR, SunlitPAR, \ - ShadedArea, SunlitArea = res.T - # 'PAR' is expected in Watt.m-2 in RATP input, whereas output is in - # micromol => convert back to W.m2 (cf shortwavebalance, line 306) - dfvox = pandas.DataFrame({'VegetationType': VegetationType, - 'Iteration': Iteration, - 'day': day, - 'hour': hour, - 'VoxelId': VoxelId, - 'ShadedPAR': ShadedPAR / 4.6, - 'SunlitPAR': SunlitPAR / 4.6, - 'ShadedArea': ShadedArea, - 'SunlitArea': SunlitArea, - 'Area': ShadedArea + SunlitArea, - 'PAR': ( - ShadedPAR * ShadedArea + SunlitPAR * - SunlitArea) / ( - ShadedArea + SunlitArea) / 4.6, - }) + dfvox = runRATP.DoIrradiation(self.ratp_grid, vegetation, sky, met) + # 'PAR' is expected in Watt.m-2 in RATP input, + # whereas output is in micromol + # => convert back to W.m2 (cf shortwavebalance, line 306) + dfvox['ShadedPAR'] /= 4.6 + dfvox['SunlitPAR'] /= 4.6 + + dfvox['Area'] = dfvox['ShadedArea'] + dfvox['SunlitArea'] + dfvox['PAR'] = (dfvox['ShadedPAR'] * dfvox['ShadedArea'] + + dfvox['SunlitPAR'] * dfvox['SunlitArea']) / dfvox['Area'] + + dfvox = dfvox.reset_index(drop=True) # for backward compatibility reasons + # TODO avoid and perform concat instead of merge return pandas.merge(dfvox, self.voxel_index()) def scene_lightmap(self, dfvox, spatial='point_id', temporal=True): diff --git a/src/alinea/pyratp/runratp.py b/src/alinea/pyratp/runratp.py index 2e0ac8a9..e48f69aa 100644 --- a/src/alinea/pyratp/runratp.py +++ b/src/alinea/pyratp/runratp.py @@ -1,147 +1,128 @@ -# Header -# - -""" -""" - -from alinea.pyratp import pyratp -import numpy as np import os +import platform import shutil import tempfile + +import numpy as np +import pandas as pd + import grid -import platform +from alinea.pyratp import pyratp + +DEBUG = False # used to debug program + +columns = ('VegetationType', + 'Iteration', + 'day', + 'hour', + 'AirTemperature', + 'VoxelId', + 'ShadedTemp', + 'SunlitTemp', + 'STARDirect', + 'STARSky', + 'ShadedPhoto', + 'SunlitPhoto', + 'ShadedTranspi', + 'SunlitTranspi', + 'ShadedArea', + 'SunlitArea', + 'ShadedGs', + 'SunlitGs', + 'ShadedAbsorbedPAR', + 'SunlitAbsorbedPAR', + 'ShadedAbsorbedNIR', + 'SunlitAbsorbedNIR') + +columns_tree = ('Iteration', + 'day', + 'hour', + 'VegetationType', + 'TotalIrradiation', + 'AirTemperature', + 'TreePhotosynthesis', + 'TreeTranspiration') + +columns_irr = ('VegetationType', + 'Iteration', + 'day', + 'hour', + 'VoxelId', + 'ShadedPAR', + 'SunlitPAR', + 'ShadedArea', + 'SunlitArea') + class runRATP(object): """ """ + @staticmethod def DoAll(*args): ratp = pyratp.ratp pyratp.dir_interception.scattering = False - ratp.out_time_spatial = np.zeros(pyratp.micrometeo.nbli*pyratp.grid3d.nveg*22*pyratp.grid3d.nent).reshape(pyratp.micrometeo.nbli*pyratp.grid3d.nveg*pyratp.grid3d.nent,22) - ratp.out_time_tree = np.zeros(pyratp.micrometeo.nbli*8*pyratp.grid3d.nent).reshape(pyratp.micrometeo.nbli*pyratp.grid3d.nent ,8) - if platform.system() == 'Windows': - path = 'c:/tmpRATP' - if os.path.exists(os.path.normpath(path)): - shutil.rmtree(os.path.normpath(path)) - os.mkdir(path) - else : - path = tempfile.mkdtemp() - os.mkdir(path+"/Resul") - grid.gridToVGX(pyratp.grid3d,path+"/Resul/","VoxelsGrid.vgx") #Save grid in VGX format - print "... grid written" - ##print np.where(pyratp.vegetation_types.ismine==1) + ratp.out_time_spatial = np.zeros( + (pyratp.micrometeo.nbli * pyratp.grid3d.nveg * pyratp.grid3d.nent, len(columns)) + ) + ratp.out_time_tree = np.zeros( + (pyratp.micrometeo.nbli * pyratp.grid3d.nent, len(columns_tree)) + ) + if DEBUG: + if platform.system() == 'Windows': + path = 'c:/tmpRATP' + if os.path.exists(os.path.normpath(path)): + shutil.rmtree(os.path.normpath(path)) + os.mkdir(path) + else: + path = tempfile.mkdtemp() + + os.mkdir(path + "/Resul") + grid.gridToVGX(pyratp.grid3d, path + "/Resul/", "VoxelsGrid.vgx") # Save grid in VGX format + print "... grid written" + + # print np.where(pyratp.vegetation_types.ismine==1) try: - numeroMin = (np.where(pyratp.vegetation_types.ismine==1))[0][0] + 1 - blMin=np.where(pyratp.grid3d.nume==numeroMin) - if len(blMin[0])>0: + numeroMin = (np.where(pyratp.vegetation_types.ismine == 1))[0][0] + 1 + blMin = np.where(pyratp.grid3d.nume == numeroMin) + if len(blMin[0]) > 0: pyratp.ratp.do_all_mine() else: pyratp.ratp.do_all() except: pyratp.ratp.do_all() - #print 'dz,', pyratp.grid3d.dz - fspatial = open(path+"/Resul"+'/spatial.txt','w') - fspatial.write('VegetationType ntime day hour AirTemperature VoxelId ShadedTemp SunlitTemp STARDirect STARSky ShadedPhoto SunlitPhoto ShadedTranspi SunlitTranspi') - fspatial.write(' ShadedArea SunlitArea ShadedGs SunlitGs ShadedAbsorbedPAR SunlitAbsorbedPAR ShadedAbsorbedNIR SunlitAbsorbedNIR') - fspatial.write('\n') - np.savetxt(fspatial,ratp.out_time_spatial,'%.6e') - fspatial.close() - ftree = open(path+"/Resul"+'/tree.txt','w') - ftree.write('ntime day hour VegetationType TotalIrradiation AirTemperature TreePhotosynthesis TreeTranspiration') - ftree.write('\n') - np.savetxt(ftree,ratp.out_time_tree,'%.6e') - ftree.close() - - # Ecriture parametres calcul_ - fichier = open(path+"/Resul"+"/data.txt", "a") - #ecriture des variable grille - fichier.write("GRILLE") - fichier.write('\n') - for d in pyratp.grid3d.__dict__: - fichier.write(str(d)) - fichier.write("\t" + str(pyratp.grid3d.__dict__[d])) - fichier.write('\n') - fichier.write('dz'+'\t') - for i in pyratp.grid3d.dz: - fichier.write(str(i) + ' ') - fichier.write('\n') - - fichier.write("VEGETATION") - fichier.write('\n') - vegetation = pyratp.vegetation_types - for jent in range(pyratp.grid3d.nent): - fichier.write('jent'+"\t") - fichier.write(str(jent)) - fichier.write('\n') - fichier.write('mu'+"\t") - fichier.write(str(vegetation.mu[jent])) - fichier.write('\n') - fichier.write('nbincli'+"\t") - fichier.write(str(vegetation.nbincli[jent])) - fichier.write('\n') - fichier.write('distinc'+"\t") - fichier.write(str(vegetation.distinc )) - fichier.write('\n') - fichier.write('nblo'+"\t") - fichier.write(str(vegetation.nblo[jent])) - fichier.write('\n') - fichier.write('nblomin'+"\t") - fichier.write(str(vegetation.nblomin)) - fichier.write('\n') - fichier.write('rf'+"\t") - fichier.write(str(vegetation.rf )) - fichier.write('\n') - fichier.write('i_gspar'+"\t") - fichier.write(str(vegetation.i_gspar[jent])) - fichier.write('\n') - fichier.write('agspar'+"\t") - fichier.write(str(vegetation.agspar[jent])) - fichier.write('\n') - fichier.write('i_gsca'+"\t") - fichier.write(str(vegetation.i_gsca[jent])) - fichier.write('\n') - fichier.write('agsca'+"\t") - fichier.write(str(vegetation.agsca[jent] )) - fichier.write('\n') - fichier.write('i_gslt'+"\t") - fichier.write(str(vegetation.i_gslt[jent])) - fichier.write('\n') - fichier.write('agslt'+"\t") - fichier.write(str(vegetation.agslt[jent])) - fichier.write('\n') - fichier.write('agsvpd'+"\t") - fichier.write(str(vegetation.agsvpd[jent])) - fichier.write('\n') - fichier.write('avcmaxn'+"\t") - fichier.write(str(vegetation.avcmaxn[jent])) - fichier.write('\n') - fichier.write('ajmaxn'+"\t") - fichier.write(str(vegetation.ajmaxn[jent])) - fichier.write('\n') - fichier.write('ardn'+"\t") - fichier.write(str(vegetation.ardn[jent])) - fichier.write('\n') - fichier.write('ismine'+"\t") - fichier.write(str(vegetation.ismine[jent])) - fichier.write('\n') - fichier.write('epm'+"\t") - fichier.write(str(vegetation.epm[jent])) - fichier.write('\n') - - - fichier.close() - - return ratp.out_time_spatial, ratp.out_time_tree + if DEBUG: + _write_result(path + "/Resul") + + # convert matrices to dataframes + df = pd.DataFrame(ratp.out_time_spatial, columns=columns) + + df['vtyp'] = [int(t) for t in df['VegetationType']] + df['iteration'] = [int(t) for t in df['Iteration']] + df['day'] = [int(t) for t in df['day']] + df['vid'] = [int(t) - 1 for t in df['VoxelId']] + + df_spatial = df.set_index(['vtyp', 'iteration', 'vid']) + + df = pd.DataFrame(ratp.out_time_tree, columns=columns_tree) + + df['vtyp'] = [int(t) for t in df['VegetationType']] + df['iteration'] = [int(t) for t in df['Iteration']] + df['day'] = [int(t) for t in df['day']] + + df_tree = df.set_index(['vtyp', 'iteration']) + + return df_spatial, df_tree @staticmethod def DoIrradiation(*args): ratp = pyratp.ratp - ratp.out_rayt = np.zeros(pyratp.micrometeo.nbli*pyratp.grid3d.nveg*9*pyratp.grid3d.nent).reshape(pyratp.micrometeo.nbli*pyratp.grid3d.nveg*pyratp.grid3d.nent ,9) + ratp.out_rayt = np.zeros( + (pyratp.micrometeo.nbli * pyratp.grid3d.nveg * pyratp.grid3d.nent, len(columns_irr)) + ) pyratp.ratp.do_interception() - + # if platform.system() == 'Windows': # path = 'c:/tmpRATP' # if os.path.exists(os.path.normpath(path)): @@ -154,4 +135,107 @@ def DoIrradiation(*args): # np.savetxt(fspatial,ratp.out_rayt,'%.6e') # fspatial.close() - return ratp.out_rayt \ No newline at end of file + # convert matrice to dataframe + df = pd.DataFrame(ratp.out_rayt, columns=columns_irr) + + df['vtyp'] = [int(t) for t in df['VegetationType']] + df['iteration'] = [int(t) for t in df['Iteration']] + df['day'] = [int(t) for t in df['day']] + df['vid'] = [int(t) - 1 for t in df['VoxelId']] + + df_rayt = df.set_index(['vtyp', 'iteration', 'vid']) + + return df_rayt + + +def _write_result(path): + # print 'dz,', pyratp.grid3d.dz + fspatial = open(path + '/spatial.txt', 'w') + fspatial.write(" ".join(columns)) + fspatial.write('\n') + np.savetxt(fspatial, pyratp.ratp.out_time_spatial, '%.6e') + fspatial.close() + + ftree = open(path + '/tree.txt', 'w') + ftree.write(" ".join(columns_tree)) + ftree.write('\n') + np.savetxt(ftree, pyratp.ratp.out_time_tree, '%.6e') + ftree.close() + + # Ecriture parametres calcul_ + fichier = open(path + "/data.txt", "a") + # ecriture des variable grille + fichier.write("GRILLE") + fichier.write('\n') + for d in pyratp.grid3d.__dict__: + fichier.write(str(d)) + fichier.write("\t" + str(pyratp.grid3d.__dict__[d])) + fichier.write('\n') + fichier.write('dz' + '\t') + for i in pyratp.grid3d.dz: + fichier.write(str(i) + ' ') + fichier.write('\n') + + fichier.write("VEGETATION") + fichier.write('\n') + vegetation = pyratp.vegetation_types + for jent in range(pyratp.grid3d.nent): + fichier.write('jent' + "\t") + fichier.write(str(jent)) + fichier.write('\n') + fichier.write('mu' + "\t") + fichier.write(str(vegetation.mu[jent])) + fichier.write('\n') + fichier.write('nbincli' + "\t") + fichier.write(str(vegetation.nbincli[jent])) + fichier.write('\n') + fichier.write('distinc' + "\t") + fichier.write(str(vegetation.distinc)) + fichier.write('\n') + fichier.write('nblo' + "\t") + fichier.write(str(vegetation.nblo[jent])) + fichier.write('\n') + fichier.write('nblomin' + "\t") + fichier.write(str(vegetation.nblomin)) + fichier.write('\n') + fichier.write('rf' + "\t") + fichier.write(str(vegetation.rf)) + fichier.write('\n') + fichier.write('i_gspar' + "\t") + fichier.write(str(vegetation.i_gspar[jent])) + fichier.write('\n') + fichier.write('agspar' + "\t") + fichier.write(str(vegetation.agspar[jent])) + fichier.write('\n') + fichier.write('i_gsca' + "\t") + fichier.write(str(vegetation.i_gsca[jent])) + fichier.write('\n') + fichier.write('agsca' + "\t") + fichier.write(str(vegetation.agsca[jent])) + fichier.write('\n') + fichier.write('i_gslt' + "\t") + fichier.write(str(vegetation.i_gslt[jent])) + fichier.write('\n') + fichier.write('agslt' + "\t") + fichier.write(str(vegetation.agslt[jent])) + fichier.write('\n') + fichier.write('agsvpd' + "\t") + fichier.write(str(vegetation.agsvpd[jent])) + fichier.write('\n') + fichier.write('avcmaxn' + "\t") + fichier.write(str(vegetation.avcmaxn[jent])) + fichier.write('\n') + fichier.write('ajmaxn' + "\t") + fichier.write(str(vegetation.ajmaxn[jent])) + fichier.write('\n') + fichier.write('ardn' + "\t") + fichier.write(str(vegetation.ardn[jent])) + fichier.write('\n') + fichier.write('ismine' + "\t") + fichier.write(str(vegetation.ismine[jent])) + fichier.write('\n') + fichier.write('epm' + "\t") + fichier.write(str(vegetation.epm[jent])) + fichier.write('\n') + + fichier.close()