Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/post_PFT.py

    r4777 r4805  
    55import datetime, netCDF4, numpy as np, os 
    66 
    7 path_luh = "/home/orchidee03/vuichard/LandUseMap_CMIP6/ESA2.0.7b_19PFT_quarterdeg_final/landusemap_%(year)s.nc" 
     7path_luh = "/home/surface7/vuichard/LandUseMap_CMIP6/ESA2.0.7b_19PFT_quarterdeg_final//landusemap_%(year)s.nc" 
    88path_orc = "/home/satellites7/vbastri/LC/orchidee/ESA2.0.7b/19PFT_quarterdeg_final/PFTmap_2010.nc" 
    9 path_out = "/home/satellites7/vbastri/LC/orchidee/ESA-LUH2-final-ext8/PFTmap_%(year)s.nc" 
     9path_out = "/home/satellites7/vbastri/LC/orchidee/ESA-LUH2-final/PFTmap_%(year)s.nc" 
    1010 
    1111# Define global attributes to be written in the output netcdf-files 
    12 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", 
    13             contact = "V.Bastrikov (vladislav.bastrikov@lsce.ipsl.fr), P.Peylin (peylin@lsce.ipsl.fr), N.Vuichard (nicolas.vuichard@lsce.ipsl.fr)", 
     12atts = 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", 
     13            contact = "V.Bastrikov (vladislav.bastrikov@lsce.ipsl.fr), P.Peylin (philippe.peylin@lsce.ipsl.fr), N.Vuichard (nicolas.vuichard@lsce.ipsl.fr)", 
    1414            web = "https://orchidas.lsce.ipsl.fr/dev/LCCCI.php", 
    1515            date = datetime.datetime.now().strftime("%d/%m/%Y %H:%M"), 
     
    2424# Load ORCHIDEE PFT data for year 2010 (with 19 PFTs) 
    2525nc = netCDF4.Dataset(path_orc) 
    26 orc19 = nc.variables["maxvegetfrac"][0] 
     26orc = nc.variables["maxvegetfrac"][0] 
    2727nc.close() 
    2828 
     
    4141    return newveg 
    4242 
    43 for year in range(1880, 2016): 
     43for year in range(850, 2016): 
    4444    print year 
    4545    nc = netCDF4.Dataset(path_luh % dict(year=year)) 
    4646    veg = nc.variables["YEARLYMAXVEGETFRAC_NORMALIZE"][0,:,::-1,:] 
    4747    nc.close 
     48 
     49    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.) 
    4850 
    4951    # Make LSM uniform across pfts 
     
    5456            print "Add %s pts to PFT%s" % (len(ys), pft+1) 
    5557            veg[pft, ys, xs] = 0 
    56     print "Total range [%s - %s], %s pts" % (veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.) 
     58    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.) 
    5759     
     60    # Clip to [0,1] and normalize 
     61    veg = np.ma.clip(veg, 0, 1) 
     62    veg = veg / veg.sum(axis=0) 
     63    print "Clip to [0,1] and normalize" 
     64    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.) 
     65 
    5866    # Fill masked inland pts with bare soil 
    59     ys, xs = np.ma.where((orc19[-2] > 0) & (veg[0].mask == True)) 
     67    ys, xs = np.ma.where((orc[-2] > 0) & (veg[0].mask == True)) 
    6068    print "Add %s pts as bare soil" % len(ys) 
    6169    veg[0, ys, xs] = 1 
    6270    for pft in range(1,15): veg[pft, ys, xs] = 0 
    63     print "Total range [%s - %s], %s pts" % (veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.) 
     71    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.) 
    6472 
    65     # Extend by 8 pixels outside the land 
    66     for idx in range(8): 
     73    # Extend by 1 pixel outside the land 
     74    for idx in range(1): 
    6775        print "Extend by 1 pixel" 
    6876        veg = extend(veg) 
    69         print "Total range [%s - %s], %s pts" % (veg.sum(axis=0).min(), veg.sum(axis=0).max(), veg.count() / 15.) 
     77        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.) 
    7078 
    7179    # Write output file 
     
    8795    ncout.variables["time_counter"].setncatts(dict(units = "years since 0-1-1", calendar = "gregorian", axis = "T")) 
    8896    ncout.variables["veget"].setncatts(dict(valid_min = 1, valid_max = len(veg), long_name = "Vegetation classes", units = "-", axis = "Z")) 
    89     ncout.variables["lat"][:] = grid.lat 
    90     ncout.variables["lon"][:] = grid.lon 
     97    ncout.variables["lat"][:] = lat 
     98    ncout.variables["lon"][:] = lon 
    9199    ncout.variables["maxvegetfrac"][0] = veg 
    92100    ncout.variables["time_counter"][0] = year 
Note: See TracChangeset for help on using the changeset viewer.