source: trunk/TOOLS/GENERATE_VEGET_MAP/ESA-LUH2/CCI_to_PFT19.py @ 6267

# $HeadURL$
# $Date$
# $Revision$

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()