Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rename modules to hydrological_processes and update imports #165

Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
66 changes: 35 additions & 31 deletions rubem/_dynamic_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,12 @@
import numpy as np
import pcraster as pcr
import pcraster.framework as pcrfw
from rubem.configuration.output_format import OutputFileFormat

from rubem.modules._evapotranspiration import *
from rubem.modules._interception import *
from rubem.modules._soil import *
from rubem.modules._surface_runoff import *
from rubem.date._date_calc import *
from rubem.file._file_generators import *
from rubem.configuration.model_configuration import ModelConfiguration
from rubem.configuration.output_format import OutputFileFormat
from rubem.date._date_calc import daysOfMonth
from rubem.file._file_generators import report
from rubem.hydrological_processes import Evapotranspiration, Interception, Soil, SurfaceRunoff


class RUBEM(pcrfw.DynamicModel):
Expand Down Expand Up @@ -77,7 +74,7 @@ def __stepReport(self):
name=outputVar,
timestep=self.currentStep,
outpath=self.config.output_directory.path,
format=OutputFileFormat.GEOTIFF,
file_format=OutputFileFormat.GEOTIFF,
base_raster_info=self.config.output_raster_base,
)

Expand Down Expand Up @@ -214,9 +211,9 @@ def initial(self):
self.ndvi_max = self.__readmap_wrapper(self.config.raster_files.ndvi_max)
self.ndvi_min = self.__readmap_wrapper(self.config.raster_files.ndvi_min)

self.logger.info("Computing min. and max. surface runoff (SR)")
self.sr_min = srCalc(self.ndvi_min)
self.sr_max = srCalc(self.ndvi_max)
self.logger.info("Computing min. and max. Reflectances Simple Ratio (SR)")
self.sr_min = Interception.get_reflectances_simple_ration(self.ndvi_min)
self.sr_max = Interception.get_reflectances_simple_ration(self.ndvi_max)

self.logger.info("Reading soil attributes...")
soil = self.__readmap_wrapper(self.config.raster_files.soil)
Expand Down Expand Up @@ -412,20 +409,20 @@ def dynamic(self):
)

self.logger.debug("Interception")
SR = srCalc(ndvi)
FPAR = fparCalc(
SR = Interception.get_reflectances_simple_ration(ndvi)
FPAR = Interception.get_fpar(
self.config.constants.fraction_photo_active_radiation_min,
self.config.constants.fraction_photo_active_radiation_max,
SR,
self.sr_min,
self.sr_max,
)
LAI = laiCalc(
LAI = Interception.get_leaf_area_index(
FPAR,
self.config.constants.fraction_photo_active_radiation_max,
self.config.constants.leaf_area_interception_max,
)
Id, Ir, Iv, self.Itp = interceptionCalc(
self.Itp = Interception.get_interception(
self.config.calibration_parameters.alpha,
LAI,
precipitation,
Expand All @@ -435,40 +432,46 @@ def dynamic(self):

self.logger.debug("Evapotranspiration")

Kc_1 = kcCalc(ndvi, self.ndvi_min, self.ndvi_max, self.kc_min, self.kc_max)
Kc_1 = Interception.get_crop_coef(
ndvi, self.ndvi_min, self.ndvi_max, self.kc_min, self.kc_max
)
# condicao do kc, se NDVI < 1.1NDVI_min, kc = kc_min
kc_cond1 = pcrfw.scalar(ndvi < 1.1 * self.ndvi_min)
kc_cond2 = pcrfw.scalar(ndvi > 1.1 * self.ndvi_min)
Kc = pcr.scalar((kc_cond2 * Kc_1) + (kc_cond1 * self.kc_min))
Ks = pcr.scalar(ksCalc(self.TUr, self.TUw, self.TUcc))
Ks = pcr.scalar(
Evapotranspiration.get_water_stress_coef_et_vegeted_area(self.TUr, self.TUw, self.TUcc)
)

# Vegetated area
self.ET_av = etavCalc(ETp, Kc, Ks)
self.ET_av = Evapotranspiration.get_et_vegetated_area(ETp, Kc, Ks)

# Impervious area
# ET impervious area = Interception of impervious area
cond = pcr.scalar((precipitation != 0))
self.ET_ai = self.config.constants.impervious_area_interception * cond

# Open water
self.ET_ao = etaoCalc(ETp, Kp, precipitation, Ao)
self.ET_ao = Evapotranspiration.get_actual_et_open_water_area(ETp, Kp, precipitation, Ao)
# self.report(ET_ao,self.config.output_directory.output+'ETao')

# Bare soil
self.ET_as = etasCalc(ETp, self.kc_min, Ks)
self.ET_as = Evapotranspiration.get_water_stress_coef_et_bare_soil_area(
ETp, self.kc_min, Ks
)
self.ETr = (Av * self.ET_av) + (Ai * self.ET_ai) + (Ao * self.ET_ao) + (As * self.ET_as)

self.logger.debug("Surface Runoff")

Pdm = precipitation / rainyDays
Ch = chCalc(
Ch = SurfaceRunoff.get_coef_soil_moist_conditions(
self.TUr,
self.dg,
self.Zr,
self.TUsat,
self.config.calibration_parameters.beta,
)
Cper = cperCalc(
Cper = SurfaceRunoff.get_runoff_coef_permeable_areas(
self.TUw,
self.dg,
self.Zr,
Expand All @@ -478,11 +481,12 @@ def dynamic(self):
self.config.calibration_parameters.w_2,
self.config.calibration_parameters.w_3,
)
Aimp, Cimp = cimpCalc(Ao, Ai)
Cwp = cwpCalc(Aimp, Cper, Cimp)
Csr = csrCalc(Cwp, Pdm, self.config.calibration_parameters.rcd)
Aimp = SurfaceRunoff.get_impervious_surface_percent_per_grid_cell(Ao, Ai)
Cimp = SurfaceRunoff.get_runoff_coef_impervious_area(Aimp)
Cwp = SurfaceRunoff.get_weighted_pot_runoff_coef(Aimp, Cper, Cimp)
Csr = SurfaceRunoff.get_actual_runoff_coef(Cwp, Pdm, self.config.calibration_parameters.rcd)

self.ES = sRunoffCalc(
self.ES = SurfaceRunoff.get_surface_runoff(
Csr,
Ch,
precipitation,
Expand All @@ -495,7 +499,7 @@ def dynamic(self):

self.logger.debug("Lateral Flow")

self.LF = lfCalc(
self.LF = Soil.get_lateral_flow(
self.config.calibration_parameters.f,
self.Kr,
self.TUr,
Expand All @@ -504,7 +508,7 @@ def dynamic(self):

self.logger.debug("Recharge Flow")

self.REC = recCalc(
self.REC = Soil.get_recharge(
self.config.calibration_parameters.f,
self.Kr,
self.TUr,
Expand All @@ -513,7 +517,7 @@ def dynamic(self):

self.logger.debug("Baseflow")

self.EB = baseflowCalc(
self.EB = Soil.get_baseflow(
self.EBprev,
self.config.calibration_parameters.alpha_gw,
self.REC,
Expand All @@ -523,7 +527,7 @@ def dynamic(self):
self.EBprev = self.EB

self.logger.debug("Soil Balance")
self.TUr = turCalc(
self.TUr = Soil.get_actual_soil_moist_cont_non_sat_zone(
self.TUrprev,
precipitation,
self.Itp,
Expand All @@ -534,7 +538,7 @@ def dynamic(self):
Ao,
self.TUsat,
)
self.TUs = tusCalc(self.TUsprev, self.REC, self.EB)
self.TUs = Soil.get_actual_water_cont_sat_zone(self.TUsprev, self.REC, self.EB)
self.TUrprev = self.TUr

self.TUsprev = self.TUs
Expand Down
8 changes: 4 additions & 4 deletions rubem/file/_file_generators.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ def report(
outpath: Union[str, bytes, os.PathLike],
base_raster_info: OutputRasterBase,
timestep: Optional[int] = None,
format: OutputFileFormat = OutputFileFormat.GEOTIFF,
file_format: OutputFileFormat = OutputFileFormat.GEOTIFF,
no_data_value: float = -9999,
):
"""Storing map data to disk using GDAL
Expand All @@ -39,16 +39,16 @@ def report(
:param name: Name used as filename. Use a filename with less than eight characters and without extension. File extension will be added automatically.
:type name: str

:param format: Output file format. Default is ``OutputFileFormat.GEOTIFF``.
:type format: OutputFileFormat, optional
:param file_format: Output file format. Default is ``OutputFileFormat.GEOTIFF``.
:type file_format: OutputFileFormat, optional

:param base_raster_info: Base raster information
:type base_raster_info: OutputRasterBase

:param no_data_value: No data value. Default is ``-9999``.
:type no_data_value: float, optional
"""
if format == OutputFileFormat.GEOTIFF:
if file_format == OutputFileFormat.GEOTIFF:
__report(
variable=variable,
timestep=timestep,
Expand Down
4 changes: 4 additions & 0 deletions rubem/hydrological_processes/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
from ._evapotranspiration import *
from ._interception import *
from ._soil import *
from ._surface_runoff import *
148 changes: 148 additions & 0 deletions rubem/hydrological_processes/_evapotranspiration.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,148 @@
import pcraster as pcr


class Evapotranspiration:
"""Class to calculate evapotranspiration.

This class contains methods to calculate potential evapotranspiration (ETp),
water stress coefficient (Ks), actual evapotranspiration of vegetated area (ETva),
pan coefficient (Kp), actual evapotranspiration of open water area (ETowa),
and actual evapotranspiration of bare soil area (ETbsa).
"""

@staticmethod
def get_water_stress_coef_et_vegeted_area(
actual_soil_moisture_content: pcr._pcraster.Field,
soil_class_wilting_point: pcr._pcraster.Field,
soild_class_field_capacity: pcr._pcraster.Field,
) -> pcr._pcraster.Field:
"""Return Water Stress Coefficient (Ks) for evapotranspiration of vegetated area.

:param actual_soil_moisture_content: Actual Soil Moisture Content [mm]
:type actual_soil_moisture_content: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param soil_class_wilting_point: Wilting Point of soil class [mm]
:type soil_class_wilting_point: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param soild_class_field_capacity: Field Capacity of soil class [mm]
:type soild_class_field_capacity: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:returns: Water Stress Coefficient (Ks) [-]
:rtype: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``
"""
ks_cond = pcr.scalar(
actual_soil_moisture_content > soil_class_wilting_point
) # when actual_soil_moisture_content < soil_class_wilting_point, (false, ks = 0)
# Multiply (actual_soil_moisture_content - soil_class_wilting_point) to avoid negative ln
return (pcr.ln((actual_soil_moisture_content - soil_class_wilting_point) * ks_cond + 1)) / (
pcr.ln(soild_class_field_capacity - soil_class_wilting_point + 1)
)

@staticmethod
def get_et_vegetated_area(
potential_etp: pcr._pcraster.Field,
crop_coef: pcr._pcraster.Field,
water_stress_coef: pcr._pcraster.Field,
) -> pcr._pcraster.Field:
"""Return evapotranspiration of vegetated area.

:param potential_etp: Potential Evapotranspiration [mm]
:type potential_etp: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param crop_coef: Crop Coefficient [-]
:type crop_coef: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param water_stress_coef: Water Stress Coefficient [-]
:type water_stress_coef: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:returns: Actual Evapotranspiration
:rtype: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``
"""
return potential_etp * crop_coef * water_stress_coef

@staticmethod
def get_pan_coef_et_open_water_area(
fetch_distance: int, wind_speed: float, relative_humidity: float
) -> pcr._pcraster.Field:
"""Return pan coefficient (Kp) for evapotranspiration of open water area.

:param fetch_distance: Fetch distance
:type fetch_distance: int

:param wind_speed: Wind speed at 2 meters [m/s-1]
:type wind_speed: float

:param relative_humidity: Relative humidity [%]
:type relative_humidity: float

:returns: pan coefficient (Kp) []
:rtype: float
"""
return (
0.482
+ 0.024 * pcr.ln(fetch_distance)
- 0.000376 * wind_speed
+ 0.0045 * relative_humidity
)

@staticmethod
def get_actual_et_open_water_area(
potential_etp: pcr._pcraster.Field,
pan_coef: pcr._pcraster.Field,
precipitation: pcr._pcraster.Field,
open_water_area_fraction: pcr._pcraster.Field,
) -> pcr._pcraster.Field:
"""Return actual evapotranspiration of open water area.

:param potential_etp: Monthly Potential Evapotranspiration [mm]
:type potential_etp: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param pan_coef: pan coefficient (Kp) []
:type pan_coef: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param precipitation: Monthly Precipitation [mm]
:type precipitation: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param open_water_area_fraction: Open water Area Fraction
:type open_water_area_fraction: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:returns: Actual evapotranspiration of open water area
:rtype: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``
"""
# condition for pixel of water
cond_water_pixel = pcr.scalar(open_water_area_fraction == 1)

partial_actual_etowa = potential_etp / pan_coef

# conditions for max value for partial_actual_etp_owa,
# if partial_actual_etp_owa > precipitation in pixel with open_water_area_fraction = 1, then actual_etp_owa = precipitation
cond_max_etp_kp = pcr.scalar((partial_actual_etowa) > precipitation)

return (
(partial_actual_etowa) * (1 - cond_water_pixel)
+ precipitation * cond_water_pixel * cond_max_etp_kp
+ (partial_actual_etowa) * cond_water_pixel * (1 - cond_max_etp_kp)
)

@staticmethod
def get_water_stress_coef_et_bare_soil_area(
potential_et: pcr._pcraster.Field,
crop_coef_min: pcr._pcraster.Field,
watet_stress_coef: pcr._pcraster.Field,
) -> pcr._pcraster.Field:
"""Return Water Stress Coefficient (Ks) for evapotranspiration of bare soil area.

:param potential_et: Monthly Potential Evapotranspiration [mm]
:type potential_et: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param crop_coef_min: Minimum crop Coefficient [-]
:type crop_coef_min: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:param watet_stress_coef: Water Stress Coefficient [-]
:type watet_stress_coef: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``

:returns: Actual Evapotranspiration of bare soil area
:rtype: pcr._pcraster.Field ``PCRASTER_VALUESCALE=VS_SCALAR``
"""
cond = 1 * pcr.scalar(watet_stress_coef != 0) # if ks is different from 0
return potential_et * crop_coef_min * watet_stress_coef * cond