Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/CCI_to_PFT19.py
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/CCI_to_PFT19.py	(revision 4776)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/CCI_to_PFT19.py	(revision 4776)
@@ -0,0 +1,223 @@
+import datetime, netCDF4, numpy as np, os
+
+path_esa   = "ESACCI-LC-L4-LCCS-Map-300m-P1Y-aggregated-0.250000Deg-2010-v2.0.7b.nc"
+path_still = "Still_et_al_2009_C4_025deg.nc"
+path_etopo = "ETOPO5.nc"
+path_out   = "PFTmap_2010.nc"
+
+# Define global attributes to be written in the output netcdf-file
+atts = dict(description = "ORCHIDEE PFT map based on the ESA CCI Land Cover map v2.0.7b aggregated by the user-tool v3.13 with the cross-walking table v2.4",
+            contact = "V.Bastrikov (vladislav.bastrikov@lsce.ipsl.fr), P.Peylin (peylin@lsce.ipsl.fr)",
+            web = "https://orchidas.lsce.ipsl.fr/dev/LCCCI.php",
+            date = datetime.datetime.now().strftime("%d/%m/%Y %H:%M"),
+            resolution = "Regular 0.25-degree",
+            pfts = "1-13: standard ORCHIDEE PFTs, 10: Temperate C3 grass, 14: Tropical C3 grass, 15: Boreal C3 grass, 16: Urban, 17: SnowIce, 18: Inland water, 19: Coastal water")
+
+# Set indices for generic PFTs
+GEN_TrBrEv  = 0
+GEN_TrBrDe  = 1
+GEN_TrNeEv  = 2
+GEN_TrNeDe  = 3
+GEN_ShBrEv  = 4
+GEN_ShBrDe  = 5
+GEN_ShNeEv  = 6
+GEN_ShNeDe  = 7
+GEN_Grass   = 8
+GEN_Crops   = 9
+GEN_BS      = 10
+GEN_Water   = 11
+GEN_SnowIce = 12
+GEN_Urban   = 13
+GEN_NoData  = 14
+GEN_LIST = ["Tree_Broadleaf_Evergreen", "Tree_Broadleaf_Deciduous", "Tree_Needleleaf_Evergreen", "Tree_Needleleaf_Deciduous", "Shrub_Broadleaf_Evergreen", "Shrub_Broadleaf_Deciduous", "Shrub_Needleleaf_Evergreen", "Shrub_Needleleaf_Deciduous", "Natural_Grass", "Crops", "Bare_soil", "Water", "Snow_Ice", "Urban", "NoData"]
+
+# Set indices for ORCHIDEE PFTs
+ORC_BS        = 0
+ORC_TropEBF   = 1
+ORC_TropDBF   = 2
+ORC_TempENF   = 3
+ORC_TempEBF   = 4
+ORC_TempDBF   = 5
+ORC_BorENF    = 6
+ORC_BorDBF    = 7
+ORC_BorDNF    = 8
+ORC_TempGrass = 9
+ORC_C4grass   = 10
+ORC_C3crops   = 11
+ORC_C4crops   = 12
+ORC_TropGrass = 13
+ORC_BorGrass  = 14
+ORC_Urban     = 15
+ORC_SnowIce   = 16
+ORC_Inland    = 17
+ORC_Coastal   = 18
+
+# Set indices for Koeppen-Geiger zones
+KG_Trop     = 1
+KG_TempWarm = 21
+KG_TempCool = 22
+KG_BorWarm  = 31
+KG_BorCool  = 32
+
+# Load ESA CCI land cover data
+nc = netCDF4.Dataset(path_esa)
+nlat = len(nc.dimensions["lat"])
+nlon = len(nc.dimensions["lon"])
+lat = nc.variables["lat"][:]
+lon = nc.variables["lon"][:]
+KG = nc.variables["user_map"][:] # Koeppen-Geiger climate zone map
+GEN = np.ma.zeros((len(GEN_LIST), nlat, nlon)) # Generic PFTs
+for idx, var in enumerate(GEN_LIST):
+    GEN[idx] = nc.variables[var][:]
+nc.close()	
+print "ESA CCI data range : [%s, %s]" % (GEN.min(), GEN.max())
+
+# Substitute masked values and/or zeros in Koeppen-Geiger map with the dominant value of each latitude or set KG_BorCool if no data
+KG = KG.filled(0)
+for ilat in range(nlat):
+    values, count = np.unique(KG[ilat], return_counts=True)
+    freq = values[np.argsort(count)][::-1]
+    if freq[0] == 0 and len(values) == 1:
+        KG[ilat] = KG_BorCool
+    else:
+        major = freq[0] if freq[0] != 0 else freq[1]
+        KG[ilat][KG[ilat]==0] = major
+print "KG data range : [%s, %s]" % (KG.min(), KG.max())
+
+# Load C4 vegetation data (Still et al., 2009)
+nc = netCDF4.Dataset(path_still)
+C4veg = nc.variables["C4veg"][:]
+nc.close()
+print "Still data range : [%s, %s]" % (C4veg.min(), C4veg.max())
+
+# Load ETOPO5 data and create land-sea mask
+etopo5 = np.zeros((2161, 4321))
+nc = netCDF4.Dataset(path_etopo)
+# Load etopo, reverse latitudes (from '-90..90' to '90..-90') and move longitudes by half (from '0..360' to '-180..180)
+etopo5[:,1:] = np.roll(nc.variables["topo"][::-1,:], 2159, axis=1)
+# Copy -180 latitude data from 180 latitude data
+etopo5[:,0] = etopo5[:,-1]
+nc.close()
+# Mask the points with height below 0
+etopo5 = np.ma.masked_where(etopo5 < 0, etopo5)
+# Create etopo mask at 0.25 resolution (if 4x4 block has at least 1 point above 0, it is considered to be land and equals to 1, otherwise equals to 0)
+etopo025 = np.zeros((nlat, nlon))
+for ilat in range(nlat):
+    for ilon in range(nlon):
+        nland = etopo5[ilat*3 : ilat*3+4, ilon*3 : ilon*3+4].count()
+        if nland > 0: etopo025[ilat, ilon] = 1
+# Add Kaspean sea
+ilat = np.where((lat >= 35) & (lat <= 50))
+ilon = np.where((lon >= 44) & (lon <= 58))
+ilon, ilat = np.meshgrid(ilon[0], ilat[0])
+etopo025[ilat, ilon] = 1
+# Add Chott lake
+ilat = np.where((lat >= 33) & (lat <= 36))
+ilon = np.where((lon >= 5) & (lon <= 7))
+ilon, ilat = np.meshgrid(ilon[0], ilat[0])
+etopo025[ilat, ilon] = 1
+print "ETOPO data at 0.25 resolution created"
+
+# Create land-sea mask with the land points where generic PFTs has anything other than water
+totveg = np.ma.sum((GEN[:11].sum(axis=0), GEN[12:].sum(axis=0)), axis=0)
+LSM = totveg.mask
+# Create inland/coastal masks
+inland_msk = np.zeros((nlat, nlon))
+coast_msk = np.zeros((nlat, nlon))
+for ilat in range(nlat):
+    for ilon in range(nlon):
+        nland = etopo025[max(0, ilat-2) : ilat+3, max(0, ilon-2) : ilon+3].sum() # calculate number of land points in 5x5 block
+        if nland == 25 and ilat < 600: # if all points in 5x5 block are land points, consider the pixel to be inland, set land-sea mask to land (not applied for Antarctica)
+            LSM[ilat, ilon] = False
+            inland_msk[ilat,ilon] = 1
+        else: # all other pixels consider to be coastal
+            coast_msk[ilat,ilon] = 1
+
+# Fill masked values in generic PFTs with zeros for easier handling
+GEN = GEN.filled(0)
+
+# Allocate Generic PFTs to ORCHIDEE PFTs
+ORC = np.zeros((19, nlat, nlon))
+
+ORC[ORC_BS] = GEN[GEN_BS]
+ORC[ORC_TropEBF] = np.where(KG==KG_Trop, GEN[GEN_TrBrEv] + GEN[GEN_ShBrEv], 0)
+ORC[ORC_TropDBF] = np.where(KG==KG_Trop, GEN[GEN_TrBrDe] + GEN[GEN_ShBrDe] + GEN[GEN_TrNeDe] + GEN[GEN_ShNeDe], 0)
+ORC[ORC_TempENF] = np.where((KG==KG_Trop)|(KG==KG_TempWarm)|(KG==KG_TempCool)|(KG==KG_BorWarm), GEN[GEN_TrNeEv] + GEN[GEN_ShNeEv], 0)
+ORC[ORC_TempEBF] = np.where((KG==KG_TempWarm)|(KG==KG_TempCool)|(KG==KG_BorWarm)|(KG==KG_BorCool), GEN[GEN_TrBrEv] + GEN[GEN_ShBrEv], 0)
+ORC[ORC_TempDBF] = np.where((KG==KG_TempWarm)|(KG==KG_TempCool)|(KG==KG_BorWarm), GEN[GEN_TrBrDe] + GEN[GEN_ShBrDe], 0) + np.where(KG==KG_TempWarm, GEN[GEN_TrNeDe] + GEN[GEN_ShNeDe], 0)
+ORC[ORC_BorENF] = np.where(KG==KG_BorCool, GEN[GEN_TrNeEv] + GEN[GEN_ShNeEv], 0)
+ORC[ORC_BorDBF] = np.where(KG==KG_BorCool, GEN[GEN_TrBrDe] + GEN[GEN_ShBrDe], 0)
+ORC[ORC_BorDNF] = np.where((KG==KG_TempCool)|(KG==KG_BorWarm)|(KG==KG_BorCool), GEN[GEN_TrNeDe] + GEN[GEN_ShNeDe], 0)
+
+# Split low vegetation to C3/C4 using Still data
+grass = GEN[GEN_Grass]
+crops = GEN[GEN_Crops]
+mix = grass + crops
+veg = GEN[GEN_TrBrEv] + GEN[GEN_TrBrDe] + GEN[GEN_TrNeEv] + GEN[GEN_TrNeDe] + GEN[GEN_ShBrEv] + GEN[GEN_ShBrDe] + GEN[GEN_ShNeEv] + GEN[GEN_ShNeDe] + GEN[GEN_Grass] + GEN[GEN_Crops]
+C4total = np.minimum(C4veg*veg, mix)
+C4grass = np.where(mix==0, 0, C4total * grass / mix)
+C3grass = grass - C4grass
+C4crops = np.where(mix==0, 0, C4total * crops / mix)
+C3crops = crops - C4crops
+
+ORC[ORC_TropGrass] = np.where(KG==KG_Trop, C3grass, 0)
+ORC[ORC_TempGrass] = np.where((KG==KG_TempWarm)|(KG==KG_TempCool), C3grass, 0)
+ORC[ORC_BorGrass]  = np.where((KG==KG_BorWarm) |(KG==KG_BorCool),  C3grass, 0)
+ORC[ORC_C4grass] = C4grass
+ORC[ORC_C3crops] = C3crops
+ORC[ORC_C4crops] = C4crops
+
+ORC[ORC_Urban] = GEN[GEN_Urban]
+ORC[ORC_SnowIce] = GEN[GEN_SnowIce]
+ORC[ORC_Inland] = GEN[GEN_Water] * inland_msk
+ORC[ORC_Coastal] = GEN[GEN_Water] * coast_msk
+
+# Limit pixels content to [0, 1] range
+print "Limit pixels out of [0,1], current min : %s, max : %s" % (ORC.min(), ORC.max())
+ORC = np.where(ORC < 0, 0, ORC)
+ORC = np.where(ORC > 1, 1, ORC)
+
+# Print statistics for ORC data
+sum_orc = ORC.sum(axis=0)
+print "ORC SUM min : %s, max : %s" % (sum_orc.min(), sum_orc.max())
+print "Number of points with sum not equal to 1 :", np.where(sum_orc != 1, 1, 0).sum()
+nout = np.where( (sum_orc > 1.0001) | (sum_orc < 0.9999), 1, 0).sum()
+print "Number of points with sum outside the interval [0.9999, 1.0001] : %s (%.4g%%)" % (nout, nout*100./(nlat*nlon))
+
+# Apply land-sea mask
+ORC = np.ma.array(ORC)
+ORC.mask = LSM
+sum_orc = ORC.sum(axis=0)
+msk = np.where(sum_orc.mask == False)
+print "ORC data after applying LSM"
+print "Number of ORC points non-masked :", sum_orc.count()
+print "Number of ORC points not equal to 1 :", np.where(sum_orc[msk] != 1, 1, 0).sum()
+nout = np.where( (sum_orc[msk] > 1.0001) | (sum_orc[msk] < 0.9999), 1, 0).sum()
+print "Number of ORC points outside the interval [0.9999, 1.0001] : %s (%.4g%%)" % (nout, nout*100/sum_orc.count())
+
+# Write file
+print "Write file", path_out
+if os.path.exists(path_out): raw_input("The file already exists! Press CTRL-C to exit or any key to continue")
+ncout = netCDF4.Dataset(path_out, "w", format = "NETCDF3_CLASSIC")
+ncout.setncatts(atts)
+ncout.createDimension("lat", nlat)
+ncout.createDimension("lon", nlon)
+ncout.createDimension("time_counter", None)
+ncout.createDimension("veget", len(ORC))
+ncout.createVariable("lat", "f8", ("lat",))
+ncout.createVariable("lon", "f8", ("lon",))
+ncout.createVariable("maxvegetfrac", "f4", ("time_counter", "veget", "lat", "lon"), fill_value = 1.e+20)
+ncout.createVariable("time_counter", "f4", ("time_counter",))
+ncout.createVariable("veget", "i4", ("veget",))
+ncout.variables["lat"].setncatts(dict(long_name="Latitude", valid_max=90., valid_min=-90., units="degrees_north", axis="Y"))
+ncout.variables["lon"].setncatts(dict(long_name="Longitude", valid_max=180., valid_min=-180., units="degrees_east", axis="X", modulo=360., topology="circular"))
+ncout.variables["maxvegetfrac"].setncatts(dict(name = "maxvegetfrac", long_name = "Vegetation types", units = "-"))
+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(ORC), long_name = "Vegetation classes", units = "-", axis = "Z"))
+ncout.variables["lat"][:] = lat
+ncout.variables["lon"][:] = lon
+ncout.variables["maxvegetfrac"][0] = ORC
+ncout.variables["time_counter"][0] = 2010
+ncout.variables["veget"][:] = range(1, len(ORC)+1)
+ncout.close()
+
Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/LC-CCI-LUT-LCCStoPFT_KG_Additional_LSCE_v2.4.csv
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/LC-CCI-LUT-LCCStoPFT_KG_Additional_LSCE_v2.4.csv	(revision 4776)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/LC-CCI-LUT-LCCStoPFT_KG_Additional_LSCE_v2.4.csv	(revision 4776)
@@ -0,0 +1,36 @@
+LCCI_Class|Koeppen_Geiger_Class|Tree_Broadleaf_Evergreen|Tree_Broadleaf_Deciduous|Tree_Needleleaf_Evergreen|Tree_Needleleaf_Deciduous|Shrub_Broadleaf_Evergreen|Shrub_Broadleaf_Deciduous|Shrub_Needleleaf_Evergreen|Shrub_Needleleaf_Deciduous_|Natural_Grass|Crops|Bare_soil|Water|Snow/Ice|Urban|NoData
+120|1|||||15|30|15||40||||||
+120|21|||||15|30|15||40||||||
+120|22|||||15|30|15||40||||||
+120|31|||||20|40|20||20||||||
+120|32|||||20|40|20||20||||||
+121|1|||||30||30||40||||||
+121|21|||||30||30||40||||||
+121|22|||||30||30||40||||||
+121|31|||||40||40||20||||||
+121|32|||||40||40||20||||||
+122|1||||||60|||40||||||
+122|21||||||60|||40||||||
+122|22||||||60|||40||||||
+122|31||||||80|||20||||||
+122|32||||||80|||20||||||
+150|1||5.0||||5.0|||20.0||70.0||||
+150|21||5.0|5.0|||5.0|5.0||30.0||50.0||||
+150|22||5.0|5.0|||5.0|5.0||30.0||50.0||||
+150|31||3.0|3.0|3.0||4.0|4.0|3.0|65.0||15.0||||
+150|32||3.0|3.0|3.0||4.0|4.0|3.0|65.0||15.0||||
+151|1||10.0|||||||20.0||70.0||||
+151|21||10.0|10.0||||||30.0||50.0||||
+151|22||10.0|10.0||||||30.0||50.0||||
+151|31||7.5|7.5|5.0|||||65.0||15.0||||
+151|32||7.5|7.5|5.0|||||65.0||15.0||||
+152|1||||||10.0|||20.0||70.0||||
+152|21||||||10.0|10.0||30.0||50.0||||
+152|22||||||10.0|10.0||30.0||50.0||||
+152|31||||||7.5|7.5|5.0|65.0||15.0||||
+152|32||||||7.5|7.5|5.0|65.0||15.0||||
+153|1|||||||||30.0||70.0||||
+153|21|||||||||50.0||50.0||||
+153|22|||||||||50.0||50.0||||
+153|31|||||||||85.0||15.0||||
+153|32|||||||||85.0||15.0||||
Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/LC-CCI-LUT-LCCStoPFT_KG_Default_LSCE_v2.4.csv
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/LC-CCI-LUT-LCCStoPFT_KG_Default_LSCE_v2.4.csv	(revision 4776)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/LC-CCI-LUT-LCCStoPFT_KG_Default_LSCE_v2.4.csv	(revision 4776)
@@ -0,0 +1,39 @@
+LCCI_Class|Tree_Broadleaf_Evergreen|Tree_Broadleaf_Deciduous|Tree_Needleleaf_Evergreen|Tree_Needleleaf_Deciduous|Shrub_Broadleaf_Evergreen|Shrub_Broadleaf_Deciduous|Shrub_Needleleaf_Evergreen|Shrub_Needleleaf_Deciduous|Natural_Grass|Crops|Bare_soil|Water|Snow/Ice|Urban|NoData
+0|||||||||||||||100
+10|||||||||10|90|||||
+11|||||||||10|90|||||
+12||||||70||||30|||||
+20|||||||||10|90|||||
+30|5|5|||5|5|5||15|60|||||
+40|7.5|7.5|||10|15|10||30|20|||||
+50|90||||5|5|||0||||||
+60||50||||20|||30||||||
+61||70||||15|||15||||||
+62||30||||25|||45||||||
+70|||70||5|5|5||15||||||
+71|||70||5|5|5||15||||||
+72|||30||||25||45||||||
+80||||50|2.5|2.5|2.5|12.5|30||||||
+81||||70|5|5|5||15||||||
+82||||30||||25|45||||||
+90||30|20|10|5|5|5||25||||||
+100|10|20|5|5|5|10|5||40||||||
+110|5|10|5||5|10|5||60||||||
+120|||||15|30|15||40||||||
+121|||||30||30||40||||||
+122||||||60|||40||||||
+130|||||||||100||||||
+140|||||||||100||||||
+150||4.7|4.7|3.1||4.7|4.7|3.1|40.0||35.0||||
+151||7.5|7.5|5.0|||||45.0||35.0||||
+152||||||7.5|7.5|5.0|45.0||35.0||||
+153|||||||||65.0||35.0||||
+160|37.5|37.5|||||||25||||||
+170|75||||25||||||||||
+180||||||25|15||60||||||
+190||||||||||||||100|
+200|||||||||||100||||
+201|||||||||||100||||
+202|||||||||||100||||
+210||||||||||||100|||
+220|||||||||||||100||
Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/post_PFT.py
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/post_PFT.py	(revision 4776)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/post_PFT.py	(revision 4776)
@@ -0,0 +1,91 @@
+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_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"
+
+# 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)",
+            web = "https://orchidas.lsce.ipsl.fr/dev/LCCCI.php",
+            date = datetime.datetime.now().strftime("%d/%m/%Y %H:%M"),
+            resolution = "Regular lat/lon 0.25-degree")
+
+# Generate latitudes/longitudes arrays and calculate areas of grid cells
+step = 0.25
+lat = np.arange(90 - step/2., -90, -step)
+lon = np.arange(-180 + step/2., 180, step)
+area = 2 * 6.371229E6**2 * np.repeat(step, len(lon))[None,:] * np.pi/180 * np.sin(0.5 * np.repeat(step, len(lat))[:,None] * np.pi/180) * np.cos(lat[:,None] * np.pi/180)
+
+# Load ORCHIDEE PFT data for year 2010 (with 19 PFTs)
+nc = netCDF4.Dataset(path_orc)
+orc19 = nc.variables["maxvegetfrac"][0]
+nc.close()
+
+# Function to extend land by 1 pixel outside land-sea mask
+def extend(veg):
+    newveg = np.ma.masked_all((len(veg), len(lat), len(lon)))
+    loc_area = np.where(veg[0].mask == False, area, 0)
+    wveg = veg * loc_area
+    for ilat in range(len(lat)):
+        for ilon in range(len(lon)):
+            if veg[0,ilat,ilon] is not np.ma.masked:
+                newveg[:,ilat,ilon] = veg[:,ilat,ilon]
+            else:
+                if wveg[:, max(0, ilat-1) : ilat+2, max(0, ilon-1) : ilon+2].count() == 0: continue
+                newveg[:,ilat,ilon] = wveg[:, max(0, ilat-1) : ilat+2, max(0, ilon-1) : ilon+2].sum(axis=-1).sum(axis=-1) / loc_area[max(0, ilat-1) : ilat+2, max(0, ilon-1) : ilon+2].sum()
+    return newveg
+
+for year in range(1880, 2016):
+    print year
+    nc = netCDF4.Dataset(path_luh % dict(year=year))
+    veg = nc.variables["YEARLYMAXVEGETFRAC_NORMALIZE"][0,:,::-1,:]
+    nc.close
+
+    # Make LSM uniform across pfts
+    vegsum = veg.sum(axis=0)
+    for pft in range(15):
+        ys, xs = np.where((vegsum.mask == False) & (veg[pft].mask == True))
+        if len(ys) > 0:
+            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.)
+    
+    # Fill masked inland pts with bare soil
+    ys, xs = np.ma.where((orc19[-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.)
+
+    # Extend by 8 pixels outside the land
+    for idx in range(8):
+        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.)
+
+    # Write output file
+    if os.path.exists(path_out % dict(year=year)): raw_input("The file already exists! Press CTRL-C to exit or any key to continue")
+    ncout = netCDF4.Dataset(path_out % dict(year=year), "w", format = "NETCDF3_CLASSIC")
+    ncout.setncatts(atts)
+    ncout.createDimension("lat", len(lat))
+    ncout.createDimension("lon", len(lon))
+    ncout.createDimension("time_counter", None)
+    ncout.createDimension("veget", len(veg))
+    ncout.createVariable("lat", "f8", ("lat",))
+    ncout.createVariable("lon", "f8", ("lon",))
+    ncout.createVariable("maxvegetfrac", "f4", ("time_counter", "veget", "lat", "lon"), fill_value = 1.e+20)
+    ncout.createVariable("time_counter", "f4", ("time_counter",))
+    ncout.createVariable("veget", "i4", ("veget",))
+    ncout.variables["lat"].setncatts(dict(long_name="Latitude", valid_max=90., valid_min=-90., units="degrees_north", axis="Y"))
+    ncout.variables["lon"].setncatts(dict(long_name="Longitude", valid_max=180., valid_min=-180., units="degrees_east", axis="X", modulo=360., topology="circular"))
+    ncout.variables["maxvegetfrac"].setncatts(dict(name = "maxvegetfrac", long_name = "Vegetation types", units = "-"))
+    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["maxvegetfrac"][0] = veg
+    ncout.variables["time_counter"][0] = year
+    ncout.variables["veget"][:] = range(1, len(veg)+1)
+    ncout.close()
+
Index: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/user-tool_run_example.txt
===================================================================
--- trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/user-tool_run_example.txt	(revision 4776)
+++ trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/user-tool_run_example.txt	(revision 4776)
@@ -0,0 +1,2 @@
+./aggregate-map.sh -PgridName=GEOGRAPHIC_LAT_LON -PnumRows=720 -PoutputLCCSClasses=true -PnumMajorityClasses=3 -PoutputPFTClasses=true -PuserPFTConversionTable=LC-CCI-LUT-LCCStoPFT_KG_Default_LSCE_v2.4.csv -PadditionalUserMap=Koeppengeiger_with_LCMap300m2010V1-4_extent_Uint8_AggAsBen.tif -PoutputUserMapClasses=true -PadditionalUserMapPFTConversionTable=LC-CCI-LUT-LCCStoPFT_KG_Additional_LSCE_v2.4.csv -PoutputAccuracy=true -PtargetDir=output ESACCI-LC-L4-LCCS-Map-300m-P1Y-2010-v2.0.7b.nc
+
