Changeset 4805 for trunk/TOOLS

Timestamp:

2017-12-01T16:10:17+01:00 (7 years ago)

Author:

josefine.ghattas

Message:

Updated to avoid negative values. Also only use 1 pixel extension. This version is used for the data post-processed on the 29 of Novembre.
Done by V Bastrikov

File:

• trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/post_PFT.py (modified) (5 diffs)

trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/post_PFT.py

@@ -5 +5 @@
 import datetime, netCDF4, numpy as np, os
 
-path_luh = "/home/orchidee03/vuichard/LandUseMap_CMIP6/ESA2.0.7b_19PFT_quarterdeg_final/landusemap_%(year)s.nc"
+path_luh = "/home/surface7/vuichard/LandUseMap_CMIP6/ESA2.0.7b_19PFT_quarterdeg_final//landusemap_%(year)s.nc"
 path_orc = "/home/satellites7/vbastri/LC/orchidee/ESA2.0.7b/19PFT_quarterdeg_final/PFTmap_2010.nc"
-path_out = "/home/satellites7/vbastri/LC/orchidee/ESA-LUH2-final-ext8/PFTmap_%(year)s.nc"
+path_out = "/home/satellites7/vbastri/LC/orchidee/ESA-LUH2-final/PFTmap_%(year)s.nc"
 
 # Define global attributes to be written in the output netcdf-files
-atts = dict(description = "ORCHIDEE PFT map based on the ESA CCI Land Cover map v2.0.7b for 2010 aggregated by the user-tool v3.13 with the cross-walking table v2.4, followed by historical reconstrution using G. Hurtt data LUH2v2h with merging rules v2.0",
-            contact = "V.Bastrikov (vladislav.bastrikov@lsce.ipsl.fr), P.Peylin (peylin@lsce.ipsl.fr), N.Vuichard (nicolas.vuichard@lsce.ipsl.fr)",
+atts = dict(description = "ORCHIDEE PFT map based on the ESA CCI Land Cover map v2.0.7b for 2010 aggregated by the user-tool v3.13 with the cross-walking table v2.4, followed by historical reconstrution using G. Hurtt data LU2v2h with merging rules v2.0",
+            contact = "V.Bastrikov (vladislav.bastrikov@lsce.ipsl.fr), P.Peylin (philippe.peylin@lsce.ipsl.fr), N.Vuichard (nicolas.vuichard@lsce.ipsl.fr)",
             web = "https://orchidas.lsce.ipsl.fr/dev/LCCCI.php",
             date = datetime.datetime.now().strftime("%d/%m/%Y %H:%M"),
@@ -24 +24 @@
 # Load ORCHIDEE PFT data for year 2010 (with 19 PFTs)
 nc = netCDF4.Dataset(path_orc)
-orc19 = nc.variables["maxvegetfrac"][0]
+orc = nc.variables["maxvegetfrac"][0]
 nc.close()
 
@@ -41 +41 @@
     return newveg
 
-for year in range(1880, 2016):
+for year in range(850, 2016):
     print year
     nc = netCDF4.Dataset(path_luh % dict(year=year))
     veg = nc.variables["YEARLYMAXVEGETFRAC_NORMALIZE"][0,:,::-1,:]
     nc.close
+
+    print "Veg range [%s - %s], total range [%s - %s], %s pts" % (veg.min(), veg.max(), veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
 
     # Make LSM uniform across pfts
@@ -54 +56 @@
             print "Add %s pts to PFT%s" % (len(ys), pft+1)
             veg[pft, ys, xs] = 0
-    print "Total range [%s - %s], %s pts" % (veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
+    print "Veg range [%s - %s], total range [%s - %s], %s pts" % (veg.min(), veg.max(), veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
     
+    # Clip to [0,1] and normalize
+    veg = np.ma.clip(veg, 0, 1)
+    veg = veg / veg.sum(axis=0)
+    print "Clip to [0,1] and normalize"
+    print "Veg range [%s - %s], total range [%s - %s], %s pts" % (veg.min(), veg.max(), veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
+
     # Fill masked inland pts with bare soil
-    ys, xs = np.ma.where((orc19[-2] > 0) & (veg[0].mask == True))
+    ys, xs = np.ma.where((orc[-2] > 0) & (veg[0].mask == True))
     print "Add %s pts as bare soil" % len(ys)
     veg[0, ys, xs] = 1
     for pft in range(1,15): veg[pft, ys, xs] = 0
-    print "Total range [%s - %s], %s pts" % (veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
+    print "Veg range [%s - %s], total range [%s - %s], %s pts" % (veg.min(), veg.max(), veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
 
-    # Extend by 8 pixels outside the land
-    for idx in range(8):
+    # Extend by 1 pixel outside the land
+    for idx in range(1):
         print "Extend by 1 pixel"
         veg = extend(veg)
-        print "Total range [%s - %s], %s pts" % (veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
+        print "Veg range [%s - %s], total range [%s - %s], %s pts" % (veg.min(), veg.max(), veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.)
 
     # Write output file
@@ -87 +95 @@
     ncout.variables["time_counter"].setncatts(dict(units = "years since 0-1-1", calendar = "gregorian", axis = "T"))
     ncout.variables["veget"].setncatts(dict(valid_min = 1, valid_max = len(veg), long_name = "Vegetation classes", units = "-", axis = "Z"))
-    ncout.variables["lat"][:] = grid.lat
-    ncout.variables["lon"][:] = grid.lon
+    ncout.variables["lat"][:] = lat
+    ncout.variables["lon"][:] = lon
     ncout.variables["maxvegetfrac"][0] = veg
     ncout.variables["time_counter"][0] = year