README: Koeppen-Geiger classfication of climate and GISS ModelE Ent Terrestrial Biosphere Model vegetation boundary conditions
AUTHOR - James Lui, Nancy Kiang
CONTACT - james.lui@nasa.gov, nancy.y.kiang@nasa.gov
To cite:  Nancy Y. Kiang, James Lui, personal communication

Scripts to automatically run and generate GISS ModelE Ent Terrestrial Biosphere Model vegetation boundary condition input files (VEG, HITEent, LAImax, LAI monthly) for alternative climates using Koeppen-Geiger (KG) classifications, following Kottek & Rubel (2006, DOI: 10.1127/0941-2948/2006/0130), with slight modification using 4-month seasons to classify As and Aw (documented in the code). 

Required modules:  R, Python, and their NetCDF libraries
See Ent_utils/README for versions.

Short summary:

(A) Lookup table:  Generate a lookup table (LUT) of Koeppen-Geiger biome classes and corresponding mixes of Ent TBM vegetation boundary conditions (lc-land cover type fraction, height-canopy height, laimax-maximum annual leaf area index, lai-monthly lai), where land cover include Ent plant fuctional types (PFTs) and bare soil cover.  This step only needs to be run once, and the user should use the frozen LUT, EntKGLUTv1.0, rather than generating their own, unless they are knowledgeable about vegetation.

(B) Vegetation boundary conditions (VBCs):  Take alternative climate and the LUT from (A) to generate alternative Ent TBM vegetation boundary conditions input files for ModelE-Ent.

(C) Paleovegetation Option:  Modify the predicted modern vegetation boundary conditions from (C) to replace with relevant cover for different time periods of the Phanerozoic.

(D) Alternative cover Option:  Take alternative map of Ent TBM land cover fractions to estimate their Koeppen-Geiger classes and associated height and LAI.  This option is suitable when an alternate Ent TBM map is obtained from a reconstruction or other vegetation model but the climate is not necessarily known.
   Rfiles/Entlc_to_KG.R:  Takes Ent TBM lc file and identifies the closest matching Koeppen-Geiger class in the lookup table generated by Step A. Generates and plots a netcdf map of the KG biomes.
        data: Ent TBM lc fractions file in format of ModelE VEG file.
        config: Example in Rfiles/config/config_Entlc_to_KG.txt, which contains a mid-Holocene reconstruction of lc, and uses the default Ent-KG LUT.

(E). R plotting scripts.  A number of R scripts can be used manually to generate maps of the various input and output files.  Check what is available in the directory Rfiles.
 
==========================================================================================================================================
Detailed summary:

A. obs2EntKGLUT.sh: 

This script generates a look up table of vegetation boundary conditions (height, LAI, LAImax, land cover fraction) based on regression of values over a Koeppen-Geiger biome. The user does not need to run this script unless they wish to generate alternative lookup tables (LUT) with other data sets. The default LUT of vegetation boundary conditions by Koeppen-Geiger biome type and Ent PFT as published is included in the repository.

  frozen LUT:  Ent_utils/kglut_modelE/EntKGLUT_HXH_fc9fd9b-2025-05-29_default.nc

This script takes observed climate and vegetation maps and generates a Koeppen-Geiger climate-biome classification of Ent PFT and bare soil pcover fraction, height, annual maximum LAI, and monthly LAI in a lookup table (LUT). Input data files and some parameter settings are specified by a configuration file. Relative paths are relative to the kglut_modelE directory.

  data: Default input observations are put in the data directory. Run data_get_HXH.sh to fetch provided data.  This script will get the data from either the NASA Center for Climate Simulation (NCCS) web data portal, or for users of the NCCS discover linux cluster, or copy from the same portal via an internal path.

  config: Default config file "config_obs2EntKGLUT_default.txt" specifies the default input observations in the data directory and parameter settings.

The config file fields are:

prec, temp: precipitation and temperature data in netCDF format for the yearly period of 12 months. The units for precipitation is mm/month and for temperature it is in Celsius. If these files have grid cells that do not have climate data or are otherwise undefined, such as some climate data that are only over land, then the biome map generated sets these grid cells as undefined . These files generate a biome file that is used by the next scripts. Note that the variable names are hard-coded to "prec" and "tmp".

lai, laimax, height, lc: data to be regressed into a lookup table by Koeppen-Geiger biome and Ent PFT. For this set, LAI must have the data for 12 months while the others are static. These files generate a lookup table in the form of csv files, one for each dataset, ordered by PFT and biome type, by regressing the data and weighing by grid surface area, and LC (for LAI, LAImax, height). It also generates another csv file which describes the number of samples taken. Some PFT/Biome combinations may not have any samples while others may only have 1. For full details on how to fine tune this part of the script, check the other README file. For LAI, the initial table may not be correct due to differences in the sampling between northern and southern hemispheres. For LC, small fractions may be present in the file.

resolution: describes the resolution of all input data. Accepted formats: 4X5, 2X2H, 1X1, HXH, QXQ. Input data must have the same resolution or it will not work.

metadata_dataversion: (optional) text describing the version of the input files.

regression: Method for regression for values of LAI, height, and LAImax. Option chooses the method of regression from weighted_average, kde (Scott's rule), or kde_silvermann (Silverman's rule). Using weighted average will calculate the mean for a given set of data using the area of the grid cell multiply by cover fraction as the weights, while KDE will find the highest density of the set of the data using the same weights (roughly equivalent to the highest peak of a histogram). This option does not affect the method of regression for land cover fractions, which will always use weighted average using the area of the grid cell as its weight..

outdir: output directory, all files generated will be placed into this directory. 

netcdf_format (NETCDF3-CLASS | NETCDF4): format of the final lookup table (see netCDF API documentation for python -> https://unidata.github.io/netcdf4-python/#creatingopeningclosing-a-netcdf-file). NETCDF3-CLASSIC is a safe option

hgt (YES|NO): YES - array names in the height dataset have the hgt_ prefix giving a netcdf name of 'hgt_<PFT>'. NO - netcdf array names do not have the hgt_ prefix.  Some datasets for modelE input have the hgt_ prepended to each PFT name for the array names.

overwrite_laimax (YES|NO): YES - LAImax field is assigned the average maximum annual LAI of the biome-PFT from the input data;  NO (default) - LAI values or remain unchanged from the input LAImax date.

lctrimfrac: threshold for land cover values to be trimmed from the final LC lookup table. Default should be 0.08.

lai_threshold_replace: threshold for monthly LAI to be replaced by the dominant hemisphere (fraction difference). Some regressed values for LAI may have a totally different range of values by hemisphere (or none at all for some colder biomes). This script aims to resolve that by having all LAI values for each hemisphere be in the same range with opposite seasonality. A value of 0 means the higher LAI values will always be used. A value of 1 means there is no replacement whatsoever. Alternatively, when lai_threshold_replace is set to any negative number, the LAI values for both hemispheres can be averaged, accounting for seasonality and relative abundance of each cover and biome. 

metadata_description: Description that should specify how the lookup table was generated (Values for cover fraction are weighed by grid surface area. Values for LAI, LAImax, and height are weighed by grid surface area and cover fraction.) as well as where they came from and the version of the files.

TO RUN:
1. Fetch input data:
   cd ${Ent_utils}/kglut_modelE/data
   ./data_get_HXH.sh
   cd ${Ent_utils}

2. Make or choose the config file
   A config file with default data is in: 
      ${Ent_utils}/kglut_modelE/config_examples/config_obs2EntKGLUT_default.txt
3. cd user
4. From the user directory, run the script with:
      ./../kglut_modelE/obs2EntKGLUT.sh <config_file>
   (NOTE: Specify the path for the <config_file>)
5. Outputs will be in the "outdir" specified in the config file.
----------------------------------------------

Calls: Rfiles/KG_classify.R, prectemp2biome.sh, regress_biome2laihite.py, trim_Ent_KG_LUT.R, csvLAIdominanthemi.py, csvremoveSTD, csv2lut.py.

KG_classify.R takes in prec and temp from the config file and generates a biome map based on the Koeppen-Geiger classification, V$(resolution)_entKG_biomes_0123abc.nc and also a pdf of the biome map, KG$(resolution)_map_XXXXXX.pdf 

regress_biome2laihite.py takes in biome from the previous script and outputs csv files with the raw data (EntKG_regression$(type)_raw_0123abc.nc) as well as a file detailing the number of samples taken (EntKG_regressionsamples_0123abc.csv)
 

Three scripts are called to filter out the csv files for the final lookup table script - trim_Ent_KG_LUT.R trims cover fraction values based on a threshold, csvLAIdominanthemi.py takes the dominant hemisphere and replaces the values of the non-dominant hemisphere in the LAI file based on a threshold of percent difference, csvremoveSTD removes the standard deviation values from the LAImax and height files.

csv2lut.py generates the final lookup table, Ent_v$(dimensions)_KoeppenGeiegerLUT_0123abc.nc, and also plots final LAI monthly lookup table

Example outputs:
KG_classify_config_<hash>_YYYY-MM-DD.txt - input config file to KG_classify.R script.
EntKGHXH_map_<hash>_YYYY-MM-DD.pdf - KG biome map
V720x360_KGbiomes_<hash>_YYYY-MM-DD.nc - netcdf file of KG biome class numbers
regress_biome2laihite_<hash>_YYYY-MM-DD.py - python script to produce regressed vegetation boundary conditions by KG biome and PFT; this is a preliminary result prior to trimming off minor PFTs found in some biomes.
EntKG_regressionLAI_monthly_raw_<hash>_YYYY-MM-DD.csv - table columns giving PFT, KG biome type, Nval|Nstd|Sval|Sstd, months 1-12, where Nval and Nstd are the NH mean and standard deviation of LAI, Sval and Sstd are the SH values.
EntKG_regressionsamples_<hash>_YYYY-MM-DD.csv - table columns giving PFT, KG biome type, number of samples, weight (cover fraction x grid area) for the hemisphere, hemisphere, where the weight is for averaging across the hemispheres.
EntKG_regressionLAI_max_raw_<hash>_YYYY-MM-DD.csv - 
EntKG_regressionheight_raw_<hash>_YYYY-MM-DD.csv
EntKG_regressionLC_raw_<hash>_YYYY-MM-DD.csv
EntKG_regressionLAI_monthly_raw_<hash>_YYYY-MM-DD.csv_LAIplot.pdf - monthly time series plots of LAI by biome and PFT.
csvLAIdominanthemi_<hash>_YYYY-MM-DD.py - python to modify the raw regresson outputs according to options in the config file for obs2EntKGLU.sh.
EntKG_regressionLAI_monthly_<hash>_YYYY-MM-DD.csv - final monthly LAI values by biome and PFT.
csvremoveSTD_<hash>_YYYY-MM-DD.py - python script to remove raw value standard deviations from the final LAImax and height tables for consistent format and clean input for generating the model input files.
EntKG_regressionheight_<hash>_YYYY-MM-DD.csv
EntKG_regressionLAI_max_<hash>_YYYY-MM-DD.csv
trim_Ent_KG_LUT_<hash>_YYYY-MM-DD.R - copy of Rfiles R script with config parameters saved to trim off small cover fractions (default < 0.08, set in config file for obs2EntKGLU.sh), remove crop cover to leave only natural vegetation, and rescale the retained PFTs to sum to 1
EntKG_regressionLC_trim_natveg_<hash>_YYYY-MM-DD.csv - table of trimmed land cover fractions.
csv2lut_<hash>_YYYY-MM-DD.py - python script to convert the final csv tables into a netcdf file.
Ent_vHXH_KoeppenGeigerLUT_<hash>_YYYY-MM-DD.nc - netcdf file of Ent Koeppen-Geiger lookup table (LUT).
Ent_vHXH_KoeppenGeigerLUT_<hash>_YYYY-MM-DD.nc_LAIplot.pdf - time series plots of LAImax and monthly LAI in the LUT by biome and PFT.

***********************************************

B. LUTaltclimate2modelEinput.sh:

This script generates input files in netCDF format based on alternate climatology to be used in modelE using input from the config file. Relative paths are relative to the kglut_modelE directory.
   LUTaltclimate2modelEinput.sh: Takes alternative climate (monthly temperature and precip) from ModelE aij files, and applies the lookup table to generate alternative ModelE map input files of vegetation boundary conditions.  
        config: User-specified, with example in "config_LUTaltclimate2modelEinput_EntKGLUTv1.0_Eocene.txt."
        calls: aij2prectemp.py, prectemp2biome.sh, lut2finalout.py, Ent_map_lc_weighted.R, Ent_maps.R

The config file fields are:

aij_dir, aij_jan: directory of the aij files to be read and the the JAN file of the aij files to be read. All 12 files JAN-DEC must exist. They should already have the correctly formatted names if output by scaleacc (XXX$(year1)-$(year2).aij$(runname).nc) and contains precipitation and temperature data. As of writing, precipitation should be in units of mm/day and temperature in Celcius.

prec, temp: (OPTIONAL) Same files as prec, temp for script obs2EntKGLUT.sh. Instead of using modelE diagnostics to generate vegetation boundary conditions, any formatted precipitation and temperature file can be used instead. If both files exist, then the step aij2prectemp.py will be skipped.

kg_biomes: (OPTIONAL) A KG biome map. Instead of using modelE diagnostics to generate vegetation boundary conditions, a KG biome map may be used instead. If the file exists, then steps aij2prectemp.py and KG_classify.R will be skipped.

lut: The path of the lookup table generated by the previous script

resolution, outdir, netcdf_format, hgt, metadata_dataversion, metadata_description: see above

TO RUN:
1. Make or choose the config file.
   An example config file with default LUT is in: 
      ${Ent_utils}/kglut_modelE/configs/config_LUTaltclimate2modelEinput_EntKGLUTv1.0_Eocene.txt
2. cd user
3. From the user directory, run the script with:
      ./../kglut_modelE/LUTaltclimate2modelEinput.sh <config_file>
   (NOTE: Specify the path of the <config_file>)
4. Outputs will be in the "outdir" specified in the config file.
----------------------------------------------

Calls:
aij2prectemp.py takes in the 12 aij files and outputs temperature and precipitation netCDF files (prec/temp_$(resolution)_$(year1)-$(year2)_$(runname)_0123abc.nc) to be read in by KG_classify.R while converting the units of precipitation from mm/day to mm/month.

KG_classify.R (see above)

lut2finalout.py takes in the biome file and the lookup table file and generates lai, laimax, height, and lc files for input into modelE based on the biome distribution and the values of the lookup table. 

Ent_map_lc_weighted.R generates lc-weighted average value maps of all pfts together for LAI, LAImax, and height and a map of the dominant pft type.

Ent_maps.R generates pft maps for all variables individually.

----------------------------------------------

Example outputs:
aij2prectemp_<hash>-YYYY-MM-DD.py - python script that is generated if user inputs ModelE aij diagnostic files for climate.
prec_4X5_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc - netcdf file of monthly climatological precipitation (mm/month) formatted for input the Rfiles R script KG_classify.R.
temp_4X5_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc - netcdf file of monthly climatological surface temperature (Celsius) formatted for input the Rfiles R script KG_classify.R.
KG_classify_config_<hash>-YYYY-MM-DD.txt - config file generated for input to KG_classify.R.
EntKG4X5_map_<hash>-YYYY-MM-DD.pdf - map of KG biome classification.
V72x46_KGbiomes_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc - netcdf file of map grid of KG biome number classes.
lut2finalout_<hash>-YYYY-MM-DD.py - python script to take the KG biome netcdf file and the LUT to generate the vegetation boundary condition netcdf files.
V72x46_lai_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc
V72x46_laimax_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc
V72x46_height_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc
V72x46_lc_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc
Ent_map_lcwtd_config_<hash>-YYYY-MM-DD.txt - config file generated for input to Rfiles script Ent_map_lc_weighted.R, which creates maps of the cominant land cover type.
V72x46_lc_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc_domlc.nc.pdf - map of dominant Ent PFT and soil land cover types.
V72x46_lc_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc_domlc.nc - netcdf file of dominant Ent PFT and soil land cover types.
V72x46_height_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD_lcwtd.pdf - map of cover-weighted grid canopy heights.
V72x46_height_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD_wtd.nc - netcdf file of cover-
V72x46_laimax_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD_lcwtd.pdf
V72x46_laimax_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD_wtd.nc
V72x46_lai_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD_lcwtd.pdf
V72x46_lai_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD_wtd.nc
V72x46_lc_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc.pdf
V72x46_height_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc.pdf
V72x46_laimax_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc.pdf
V72x46_lai_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD.nc.pdf
V72x46_lc_3507-3516_P2SAoM40_Phan85Ma1500_2_<hash>-YYYY-MM-DD_85.0Ma.nc

***********************************************

B-1. modelE2KGBiomes.sh

This script makes it easier to convert modelE diagnostics files into a Koeppen-Geiger map. This is essentially just the first two steps of the LUTaltclimate2modelEinput.sh script. Its input config file, config_modelE2KGBiomes.txt is the same config file as config_LUTaltclimate2modelEinput

***********************************************
B-2. paleo_entcover.py

This script takes any set of vegetation boundary conditions netcdf files  (cover, LAI maximum, height, monthly LAI) and replaces modern vegetation with the relevant paleo-vegetation, given the input year in Mya.  It runs from the command line.  Below is a template command, in which the file names can be changed to the user's file names.

python ../../kglut_modelE/paleo_entcover.py  -i V72x46_lc_${RUN}_${DATE}.nc -o V72x46_lc_${RUN}_paleo.nc -ilm V72x46_laimax_${RUN}_${DATE}.nc -olm V72x46_laimax_${RU
N}_paleo.nc -il V72x46_lai_${RUN}_${DATE}.nc -ol V72x46_lai_${RUN}_paleo.nc -ih V72x46_height_${RUN}_${DATE}.nc -oh V72x46_height_${RUN}_paleo.nc -t {$Mya}

***********************************************

C. Rfiles/Entlc_to_KG.R

This R script takes an Ent TBM lc file and identifies the closest matching Koeppen-Geiger class in the lookup table generated by Step A. It uses simple minimization of sums of square differences between existing lc fractions and the mix of Ent TBM lc fractions in the KG LUT.  Generates and plots a netcdf map of the KG biomes.
        data: Ent TBM lc fractions file in format of ModelE VEG file.
        config: Example in Rfiles/config/config_Entlc_to_KG.txt, which contains a mid-Holocene reconstruction of lc, and uses the default Ent-KG LUT.

***********************************************

D. R plotting scripts.  It is best to scan through the R scripts in the directory Rfiles to see what is available.  These generaly take input config files, for which examples are provided in Rfiles/config.

To see instructions how to use each script:
   Rscript <R script file name>


