diff --git a/couplage/MAR_LU2CARAIB.py b/couplage/MAR_LU2CARAIB.py new file mode 100644 index 0000000000000000000000000000000000000000..be3b5ef42eca0c4347de16979e28830c854f2424 --- /dev/null +++ b/couplage/MAR_LU2CARAIB.py @@ -0,0 +1,131 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Jul 10 14:44:05 2023 + +@author: Thomas +""" + +# +# ---------------------------------------------------------------------------------------------------------------------- +# +# Title of script +# +# ---------------------------------------------------------------------------------------------------------------------- + + + +# ---------------------------------------------------------------------------------------------------------------------- +# +# Imports and global variables +# \**********************************/ +# +import numpy as np + +import time + +import xarray as xr + +import os + +import sys +import matplotlib.pyplot as plt +# ---------------------------------------------------------------------------------------------------------------------- +# +# Functions +# \***************/ +# +# +# Here you can define usefull functions to be used in the main +# +"""##############################################################################################################""" +def convert(x): + converteur = {0:5, # Barren soil + 1:1, # Crops low + 2:1, # Crops medium + 3:1, # Crops high + 4:2, # Grass low + 5:2, # Grass medium + 6:2, # Grass high + 7:0, # Broadleaf low + 8:0, # Broadleaf medium + 9:0, # Broadleaf high + 10:0, # Needleleaf low + 11:0, # Needleleaf medium + 12:0, # Needleleaf high + 13:3, # City + } + + + return converteur[x] + + + + +# Main +# \**********/ +# +# +# Here you can define the instructions that are called when you execute this file +# +if __name__ == '__main__': + t0=time.time() + + # Convert degrees to radians for the central longitude and latitude + + + mar = xr.open_dataset("E:/MAR_data_raw_BE/MERRA/ICE_a01_2012.nc") + + mar = mar.isel(time=1) + + longitude = mar.LON.values + latitude = mar.LAT.values + + veg = mar.VEG.values + + frac = mar.FRV.values + + sect = [0,1,2] + + f = open('./test/land_use_caraib.dat', "w") + + print(range(len(mar.x.values))) + for i in range(len(mar.x.values)): + for j in range(len(mar.y.values)): + + veg_car = np.zeros(6) + + lo = longitude[j,i] + la = latitude[j,i] + + new_lu = [convert(veg[x,j,i]) for x in sect] + + new_frv = [frac[x,j,i]/100 for x in sect] + + if not np.isclose(np.sum(new_frv),1,atol=0.0001): + print('error') + print(np.sum(new_frv)) + new_frv = [(frac[x,j,i]/100)/np.sum(frac[x,j,i]/100) for x in sect] + print(new_frv) + sys.exit() + + for k,values in enumerate(new_lu): + + veg_car[values]= veg_car[values] + new_frv[k] + + + + veg_car = np.round(veg_car,2) + + if not np.isclose(np.sum(veg_car),1,atol=0.0001): + difference = 1 - np.sum(veg_car) + max_index = np.argmax(veg_car) + veg_car[max_index] += difference + + # print(lo,la) + # sys.exit() + text ='\t' + '{:.4f}'.format(round(lo,4))+'\t' + '{:.4f}'.format(round(la,4)) +'\t'+ '\t'.join('{:.2f}'.format(round(x,2)) for x in veg_car) + "\n" + f.write(text) + + f.close() + t1=time.time() + print('Programme compiled in ',t1-t0, 'sec') \ No newline at end of file