Upload New File

couplage/MAR_LU2CARAIB.py

+# -*- 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