diff --git a/README.md b/README.md index 1cca1be..82ac768 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,48 @@ -# GROWTH SIMULATOR +# STAND GROWTH SIMULATOR +### SIMULADOR DE CRECIMIENTO EN RODALES -A timber plantation growth simulator, with stands and management policies options. +This repo contains a timber plantation growth simulator, for stands generating various management policy options. It's entirely made in python with simple dependencies. -Based on !["Modelos de predicción de biomasa a nivel de rodal en plantaciones de Eucalyptus globulus y Pinus radiata en Zona centro sur de en Chile"](Modelos%20de%20predicción%20de%20biomasa%20a%20nivel%20de%20rodal%20en%20plantaciones%20de%20Eucalyptus%20globulus%20en%20Chile.pdf). By: Alejandro Miranda, Blas Mola and Víctor Hinojosa +The code is based on the paper: !["Modelos de predicción de biomasa a nivel de rodal en plantaciones de Eucalyptus globulus y Pinus radiata en Zona centro sur de en Chile"](Modelos%20de%20predicción%20de%20biomasa%20a%20nivel%20de%20rodal%20en%20plantaciones%20de%20Eucalyptus%20globulus%20en%20Chile.pdf). By: Alejandro Miranda, Blas Mola and Víctor Hinojosa -1. For 34 types of eucalyptus and pinus plantations in the central-south of Chile, a statistical study fitted the following power law (tabla.csv): - +For 34 types of eucalyptus and pinus plantations in the central-south of Chile, a statistical study fitted the power law found in `tabla.csv`. +The result file of the statistical study includes the following information for each of the 34 combinations of species and rotation type: + +|Attribute|Translation description| +|--|--| +|id| model identification representing the type of stand| +|Next| id of the model the stand would turn into if it is thinned. When `null`, the stand can't be thinned (eucalyptus for example).| +|Prev|Previous, opposite of next; the id of the stand before it would've been thinned.| +|Especie|species of the stand (pinus or eucalyptus).| +|Zona de crecimiento|Zone of Chile where the stand is growing. It depends on the species of the stand (pinus and eucalyptus have diferent zones).| +|Edad del rodal|Age of the stand.| +|Densidad Inicial (n init)|Initial density of the stand for ha.| +|Índice de Sitio (SI)|Site index, mean height of the stands of 10 years.| +|Manejo aplicado al rodal|management of the stand (Pulpable, Multiproposito, Intensivo 1, Intensivo 2 or Sin manejo).| +|Condición|Condition of the stand. It can be "posterior Raleo" (post thinning), "posterior a poda y raleo" (post thinning and prunning)| +|$\alpha$,$\beta$, $\gamma$|coeficents for the biomass model| +|stableyear|year when the formula begins to yield stable results (it depends on the type of pinus).| + + + +The managment of the stand, can be anyone of the next: + +• ***Pulpable***: pulpwood, case where the stand thinn for disposal. + +• ***Multipropósito***: multiporpuse, when the stand is thinned and pruned 2 times. + +• ***Intensivo 1***: Intensive 1, when the stand pruned 3 times and thinned 2 times. + +• ***Intensivo 2***: Intensive 2, like intensive 1, but with a different site cuality. + +• ***Sin manejo***: wihtout any managment.| + +1. The equation to get biomass for a stand of age t is: $$ biomass(t) = \alpha \cdot t^\beta + \gamma $$ -1. To extrapolate to earlier years; due to the formula yielding negative values (for a few types of pinus). A linear adjustment was made: +2. To extrapolate to earlier years; due to the formula yielding negative values (for a few types of pinus). A linear adjustment was made: $$ \text{biomass}(t) = @@ -22,10 +54,11 @@ $$ _Stable year is the year when the formula begins to yield stable results (it depends on the type of pinus)._ -1. A template for generating a timber plantation and different management policies was made (config.toml) +3. A template for generating a timber plantation and different management policies was made (config.toml) + ```toml - horizonte = 10 # number of years to generate - rodales = 36 # number of stands to generate, choosing one model at random + horizon = 10 # number of years to generate + [random] # random number generator seed: omit or comment for different results each run @@ -33,21 +66,22 @@ _Stable year is the year when the formula begins to yield stable results (it dep # `low` (inclusive) to `high` (exclusive) # n, n+1 for getting a single value: n # see np.random.randint - edades = [1, 18] # min, max age of generated stands + ages = [1, 18] # min, max age of generated stands has = [5, 15] # min, max hectares of generated stands + stands = 36 # number of stands to generate, choosing one model at random # ranges: start, stop, step # n, n+1 for getting a single value range: [n] # see np.arange [eucalyptus] - cosechas = [10, 20, 3] # for each Eucalyptus stand, generate different biomass histories harvesting in the year 10, 13, 16 and 19 (4 histories) + harvest = [10, 20, 3] # for each Eucalyptus stand, generate different biomass histories harvesting in the year 10, 13, 16 and 19 (4 histories) - [pino] - raleos = [6, 11, 2] # thinning policies in the years 6, 8, 10. - cosechas = [18, 29, 3] # harvesting policies in the years 18, 21, 24, ... (every 3 years) + [pino] #pinus + thinning = [6, 11, 2] # thinning policies in the years 6, 8, 10. + harvest = [18, 29, 3] # harvesting policies in the years 18, 21, 24, ... (every 3 years) # all feasible histories combining thinning and harvesting policies will be generated - min_ral = 6 #lower bound to tree thinning actions - ``` + min_thinning = 6 #lower bound to tree thinning actions + ``` 1. For a real instance, the numbers of stands and its characteristics can be passed instead of generated at random. For this a `.csv` or `.shp` file is used, the data must include the following fields: @@ -56,7 +90,7 @@ _Stable year is the year when the formula begins to yield stable results (it dep - **age**: Age of the stand - **hectare (ha)**: Area in hectares -See `auxiliary.py` for GIS extraction or creation of this attributes table `bosque_data.csv` using the methods `get_data` and `create_bosque`. +See `auxiliary.py` for GIS extraction or creation of this attributes table `bosque_data.csv` using the methods `get_data` and `create_forest`. ### quick start diff --git a/auxiliary.py b/auxiliary.py index 680239a..80e22d0 100644 --- a/auxiliary.py +++ b/auxiliary.py @@ -48,11 +48,12 @@ def get_data(filepath=".\\test\\data_base\\proto.shp"): datasource = None""" -def create_forest(gdf, id="fid"): - """Crea un bosque a partir de un geopandas dataframe""" +def create_forest(gdf, id="fid", mid="growth_mid", outfile="bosque_data.csv"): + if "area_ha" not in gdf.columns: + gdf["area_ha"] = gdf.geometry.area / 1e4 data_rodales = gdf.dropna(subset=["edad"]) data_rodales_2 = data_rodales.loc[data_rodales["area_ha"] > 0] - bos_names = ["rid", "mid", "edad_inicial", "ha"] # aprender hacer formato decente + bos_names = ["rid", "growth_mid", "edad_inicial", "ha"] # aprender hacer formato decente rodales = [] for idx, r in data_rodales_2.iterrows(): @@ -61,8 +62,8 @@ def create_forest(gdf, id="fid"): e0 = r["edad"] ha = r["area_ha"] rodal = { - "rid": r[id], # r["fid"] en caso habitual - "mid": r["id"], + "rid": r[id], + "growth_mid": r[mid], "edad_inicial": e0, "ha": ha, } @@ -70,11 +71,9 @@ def create_forest(gdf, id="fid"): bos = np.array( [tuple(r[k] for k in bos_names) for r in rodales], - dtype=[("rid", "i4"), ("mid", "i4"), ("edad_inicial", "i4"), ("ha", "f4")], - ) - np.savetxt( - "bosque_data.csv", bos, delimiter=",", header=",".join(bos_names), comments="", fmt=["%d", "%d", "%d", "%.2f"] + dtype=[("rid", "i4"), ("growth_mid", "i4"), ("edad_inicial", "i4"), ("ha", "f4")], ) + np.savetxt(outfile, bos, delimiter=",", header=",".join(bos_names), comments="", fmt=["%d", "%d", "%d", "%.2f"]) def plot_1_id_model(horizon: int = 40, show=True, save=False, target_id: int = 30): diff --git a/config.toml b/config.toml index 3569e3a..54b2f19 100644 --- a/config.toml +++ b/config.toml @@ -1,5 +1,11 @@ -horizonte = 10 -rodales = 36 +# +horizon = 10 +#tabla de los modelos de crecieminto +models_table="tabla.csv" +study_area = "fdo/instance1_growth_attributes.shp" +id = "fid" + + [random] # random number generator seed: omit or comment for different results each run @@ -7,19 +13,18 @@ seed = 4 # `low` (inclusive) to `high` (exclusive) # n, n+1 for getting a single value: n # see np.random.randint -edades = [1, 18] +ages = [1, 18] has = [5, 15] +stands = 36 # ranges: start, stop, step # n, n+1 for getting a single value range: [n] # see np.arange [pino] -podas = [5, 15, 5] # not implemented -raleos = [6, 11, 1] -cosechas = [18, 29, 1] -min_ral = 6 +pruning = [5, 15, 5] # not implemented +thinning = [6, 11, 1] +harvest = [18, 29, 1] +min_thinning = 6 [eucalyptus] -podas = [4, 6, 2] # not implemented -cosechas = [8, 13, 1] - +harvest = [8, 13, 1] diff --git a/simulator.py b/simulator.py index a7e2cbb..fbc960f 100644 --- a/simulator.py +++ b/simulator.py @@ -6,13 +6,13 @@ rodales = [ ... {'rid': 9, # rodal id - 'mid': 24, # model id + 'growth_mid': 24, # model id 'edad_inicial': 17, 'edad_final': 27, 'ha': 14, 'manejos' : [ ... {'rid': 9, - 'cosecha': 18, + 'harvest': 18, 'raleo': 6, 'biomass': array([2593., 0., 54.7, 76.9, 118.8, 182.7, 270.2, 40., 46.7, 53.4]), 'edades': array([17, 0, 1, 2, 3, 4, 5, 6, 7, 8]), @@ -41,7 +41,7 @@ - get_models: leer modelos de crecimiento desde un archivo csv - read_toml: leer configuración desde un archivo toml - calc_biomass: calcular la biomasa para un model y una edad - - generar_codigo_kitral: generar un diccionario de códigos Kitral basado en la Especie, edad y condición + - generate_kitral_code: generar un diccionario de códigos Kitral basado en la Especie, edad y condición - write: escribir archivos de salida - print_manejos_possibles: imprimir los manejos posibles @@ -53,6 +53,7 @@ Notice: Numpy is set to print only one decimal digit """ + import sys from itertools import product from pathlib import Path @@ -105,7 +106,7 @@ def calc_biomass(model: np.void, e: int) -> float: return model["α"] * e ** model["β"] + model["γ"] -def generar_codigo_kitral(especie: str, edad: int, condicion: str) -> int: +def generate_kitral_code(especie: str, edad: int, condicion: str) -> int: """Genera un diccionario de códigos Kitral basado en la Especie, edad y condición""" if especie == "pino": if edad <= 3: @@ -129,23 +130,65 @@ def generar_codigo_kitral(especie: str, edad: int, condicion: str) -> int: return value +def get_kitral_data(codigo_kitral: int) -> dict: + """Obtiene la especie, rango de edad y condición a partir de un código Kitral.""" + if codigo_kitral in [19, 20, 21, 22, 23, 24, 25]: + especie = "pino" + if codigo_kitral == 19: + edad = [0, 4] + condicion = "PostRaleo1250-700" + elif codigo_kitral == 20: + edad = [4, 12] + condicion = "PostRaleo1250-700" + elif codigo_kitral == 21: + edad = [12, 18] + condicion = "PostRaleo1250-700" + elif codigo_kitral == 22: + edad = [18, 41] + condicion = "PostRaleo1250-700" + elif codigo_kitral == 23: + edad = [6, 12] + condicion = "PostPodayRaleo700-300" + elif codigo_kitral == 24: + edad = [12, 18] + condicion = "PostPodayRaleo700-300" + elif codigo_kitral == 25: + edad = [18, 41] + condicion = "PostPodayRaleo700-300" + elif codigo_kitral in [26, 27, 28]: + especie = "eucalyptus" + if codigo_kitral == 26: + edad = [0, 4] + condicion = "SinManejo" + elif codigo_kitral == 27: + edad = [4, 11] + condicion = "SinManejo" + elif codigo_kitral == 28: + edad = [11, 31] + condicion = "SinManejo" + else: + return {"error": "Código Kitral desconocido"} + + return {"especie": especie, "edad": edad, "condicion": condicion} + + def print_manejos_possibles(config): """Imprime todos los manejos posibles para los rodales""" manejos_posibles = [] print("manejos posibles", end=": ") - for cosecha, raleo in product(np.arange(*config["pino"]["cosechas"]), np.arange(*config["pino"]["raleos"])): + for cosecha, raleo in product(np.arange(*config["pino"]["harvest"]), np.arange(*config["pino"]["thinning"])): if raleo > cosecha: continue print(f"(c{cosecha}, r{raleo})", end=", ") manejos_posibles.append([int(raleo), int(cosecha)]) - for cosecha in np.arange(*config["eucalyptus"]["cosechas"]): + for cosecha in np.arange(*config["eucalyptus"]["harvest"]): print(f"(c{int(cosecha)}, r{-1})", end=", ") manejos_posibles.append([-1, int(cosecha)]) - for cosecha in np.arange(*config["pino"]["cosechas"]): + for cosecha in np.arange(*config["pino"]["harvest"]): print(f"(c{cosecha}, r{-1})", end=", ") manejos_posibles.append([-1, int(cosecha)]) - for raleo in np.arange(*config["pino"]["raleos"]): + for raleo in np.arange(*config["pino"]["thinning"]): print(f"(c{-1}, r{raleo})", end=", ") manejos_posibles.append([int(raleo), -1]) @@ -166,8 +209,7 @@ def read_toml(config_toml="config.toml"): return config -def generate_random_forest(config=read_toml(), models=get_models()): - +def generate_random_forest(config=None, models=None): # 0 setup random number generator if seed := config["random"].get("seed"): rng = np.random.default_rng(seed) @@ -178,16 +220,16 @@ def generate_random_forest(config=read_toml(), models=get_models()): rodales = [] # itera = iter(range(config["rodales"])) # r = next(itera) - for r in range(config["rodales"]): + for r in range(config["random"]["stands"]): model = rng.choice(models) # model = rng.choice(models) # print(model) - e0 = rng.integers(*config["random"]["edades"]) - e1 = e0 + config["horizonte"] + e0 = rng.integers(*config["random"]["ages"]) + e1 = e0 + config["horizon"] ha = rng.integers(*config["random"]["has"]) rodal = { - "rid": rodal["rid"], - "mid": model["id"], + "rid": r, + "growth_mid": model["id"], "edad_inicial": e0, "edad_final": e1, "ha": ha, @@ -197,17 +239,16 @@ def generate_random_forest(config=read_toml(), models=get_models()): return rodales -def generate_forest(config=read_toml(), filepath="./bosque_data.csv"): - +def generate_forest(config=None, filepath="bosque_data.csv"): data = np.genfromtxt(filepath, delimiter=",", names=True) rodales = [] for r in data: e0 = r["edad_inicial"] - e1 = e0 + config["horizonte"] + e1 = e0 + config["horizon"] ha = r["ha"] rodal = { "rid": r["rid"], # Identificador único para cada rodal - "mid": r["mid"], + "growth_mid": r["growth_mid"], "edad_inicial": e0, "edad_final": e1, "ha": ha, @@ -217,10 +258,10 @@ def generate_forest(config=read_toml(), filepath="./bosque_data.csv"): return rodales -def generate(config=read_toml(), models=get_models(), rodales=generate_forest()): +def generate(config=None, models=None, rodales=None): """Genera los rodales con las biomasas generadas por cada año, dependiendo de su manejo y edad de crecimiento, junto con la biomasa para vender y el codigo kitral""" for rodal in rodales: - indices = np.where(models["id"] == rodal["mid"])[0] + indices = np.where(models["id"] == rodal["growth_mid"])[0] model = models[indices][0] # model = rng.choice(models) # print(model) @@ -237,33 +278,33 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) "edades": edades, "eventos": ["" for e in edades], "vendible": [0 for e in edades], - "codigo_kitral": [generar_codigo_kitral(model["Especie"], e, "sin manejo") for e in edades], + "codigo_kitral": [generate_kitral_code(model["Especie"], e, "sin manejo") for e in edades], } ] - # has cosecha if any of the proposed "cosechas" ranges are in the simulated "edades" - has_cosecha = any(np.isin(np.arange(*config[model["Especie"]]["cosechas"]), edades)) + # has cosecha if any of the proposed "harvest" ranges are in the simulated "edades" + has_cosecha = any(np.isin(np.arange(*config[model["Especie"]]["harvest"]), edades)) # can have raleo only if it's pino if model["Especie"] == "pino": if not has_cosecha: # easy case - has_raleo = (model["next"] != -1) and any(np.isin(np.arange(*config["pino"]["raleos"]), edades)) + has_raleo = (model["next"] != -1) and any(np.isin(np.arange(*config["pino"]["thinning"]), edades)) else: # adjust "edades" -> "edades_manejo", by periodically "cosecha" (to harvest) via modulus operator # then check if raleo is in range has_raleo = False - for cosecha in np.arange(*config["pino"]["cosechas"]): + for cosecha in np.arange(*config["pino"]["harvest"]): edades_manejo = edades % cosecha if model["prev"] == -1: # no raleado desde un inicio has_raleo = (model["next"] != -1) and any( - np.isin(np.arange(*config["pino"]["raleos"]), edades_anejo) + np.isin(np.arange(*config["pino"]["thinning"]), edades_manejo) ) else: # raleado desde un inicio - raleos = np.arange(*config["pino"]["raleos"]) - has_raleo = any((raleo + cosecha) in edades for raleo in raleos) + thinning = np.arange(*config["pino"]["thinning"]) + has_raleo = any((raleo + cosecha) in edades for raleo in thinning) - print(cosecha, np.arange(*config["pino"]["raleos"]), edades_manejo, has_raleo) + print(cosecha, np.arange(*config["pino"]["thinning"]), edades_manejo, has_raleo) if has_raleo: break else: @@ -277,9 +318,9 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) display(manejos) # 2 elif has_cosecha and not has_raleo: - # iterb = iter(np.arange(*config[model["Especie"]]["cosechas"])) + # iterb = iter(np.arange(*config[model["Especie"]]["harvest"])) # cosecha = next(iterb) - for cosecha in np.arange(*config[model["Especie"]]["cosechas"]): + for cosecha in np.arange(*config[model["Especie"]]["harvest"]): if cosecha not in edades: display(f"skipping: {e0=} !< {cosecha=} !< {e1=}") continue @@ -288,7 +329,8 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) mods = [model["id"]] * len(edades_manejo) else: mods = [ - model["id"] if e > config[model["Especie"]]["min_ral"] else model["prev"] for e in edades_manejo + model["id"] if e > config[model["Especie"]]["min_thinning"] else model["prev"] + for e in edades_manejo ] manejo = { "rid": rodal["rid"], @@ -300,9 +342,9 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) "vendible": ha * np.array([calc_biomass(model, cosecha) if e == 0 else 0 for e in edades_manejo]), "codigo_kitral": [ ( - generar_codigo_kitral(model["Especie"], cosecha, "sin manejo") + generate_kitral_code(model["Especie"], cosecha, "sin manejo") if e == 0 - else generar_codigo_kitral(model["Especie"], e, "sin manejo") + else generate_kitral_code(model["Especie"], e, "sin manejo") ) for e in edades_manejo ], @@ -311,9 +353,9 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) display(manejo) # 3 elif not has_cosecha and has_raleo: - # iterc = iter(np.arange(*config[model["Especie"]]["raleos"])) + # iterc = iter(np.arange(*config[model["Especie"]]["thinning"])) # raleo = next(iterc) - for raleo in np.arange(*config[model["Especie"]]["raleos"]): + for raleo in np.arange(*config[model["Especie"]]["thinning"]): if raleo not in edades: display(f"skipping: {e0=} !< {raleo=} !< {e1=}") continue @@ -357,9 +399,9 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) ), "codigo_kitral": [ ( - generar_codigo_kitral(model["Especie"], e, "sin manejo") + generate_kitral_code(model["Especie"], e, "sin manejo") if e < raleo - else generar_codigo_kitral(model["Especie"], e, "con manejo") + else generate_kitral_code(model["Especie"], e, "con manejo") ) for e in edades ], @@ -370,16 +412,15 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) else: # has cosecha and raleo, se asume que se raleo antes del periodo 0 en calc_biomass # iterd = iter( # product( - # np.arange(*config[model["Especie"]]["cosechas"]), np.arange(*config[model["Especie"]]["raleos"]) + # np.arange(*config[model["Especie"]]["harvest"]), np.arange(*config[model["Especie"]]["thinning"]) # ) # ) # cosecha, raleo = next(iterd) for cosecha, raleo in product( - np.arange(*config[model["Especie"]]["cosechas"]), np.arange(*config[model["Especie"]]["raleos"]) + np.arange(*config[model["Especie"]]["harvest"]), np.arange(*config[model["Especie"]]["thinning"]) ): edades_manejo = edades % cosecha if model["prev"] == -1: - if (raleo >= cosecha) or (cosecha not in edades) or (raleo not in edades_manejo): # display(f"skipping: {min(edades_manejo)=} {max(edades_manejo)=} !< {raleo=} !< {cosecha=} !< {e1=}") continue @@ -462,12 +503,12 @@ def generate(config=read_toml(), models=get_models(), rodales=generate_forest()) "vendible": ha * np.array(vendible), "codigo_kitral": [ ( - generar_codigo_kitral(model["Especie"], e, "con manejo") + generate_kitral_code(model["Especie"], e, "con manejo") if e >= raleo else ( - generar_codigo_kitral(model["Especie"], cosecha, "con manejo") + generate_kitral_code(model["Especie"], cosecha, "con manejo") if e == 0 - else generar_codigo_kitral(model["Especie"], e, "sin manejo") + else generate_kitral_code(model["Especie"], e, "sin manejo") ) ) for e in edades_manejo @@ -496,7 +537,7 @@ def write(rodales): np.savetxt("vendible.csv", vd.T, delimiter=",", header=names, comments="") np.savetxt("codigo_kitral.csv", kt.T, delimiter=",", header=names, comments="", fmt="%s") - bos_names = ["rid", "mid", "edad_inicial", "ha"] # aprender hacer formato decente + bos_names = ["rid", "growth_mid", "edad_inicial", "ha"] # aprender hacer formato decente bos = np.array([tuple(r[k] for k in bos_names) for r in rodales]) np.savetxt("bosque.csv", bos, delimiter=",", header=",".join(bos_names), comments="", fmt="%d") @@ -577,7 +618,6 @@ def main(argv=None): rodales_sin_manejo = generate_random_forest() else: - # usar bosque_data.csv, si no se tiene se puede crear con las funciones del auxiliary rodales_sin_manejo = generate_forest(config, args.data_forest)