aquacropos / aquacrop

AquaCrop-OSPy: Python implementation of AquaCrop-OS
https://aquacropos.github.io/aquacrop/
Apache License 2.0
95 stars 68 forks source link

biomass_accumulation doesn't handle Et0=0.0 #2

Closed wojeda-mx closed 3 years ago

wojeda-mx commented 3 years ago

When running the brussels_climate.txt included dataset, there is a datapoint with ET=0 on 16/12/1980: 16 12 1980 1.7 6.9 0.2 0.0 .

This raises a ZeroDivisionError when running model.step on biomass_accumulation. We believe it's because of the Et0 denominator on https://github.com/thomasdkelly/aquacrop/blob/f4c489dcb5073096047e09d60aac7658ec1f2bd0/aquacrop/solution.py#L3732-L3734 , and we could either try/catch, or use numpy to cast the floats and ensure we get NaN/Inf for such cases (https://stackoverflow.com/questions/10011707/how-to-get-nan-when-i-divide-by-zero).

tested using the collab Notebook 1, only switching the dataset. stack trace from %xmode Verbose

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-113-6829cdd20957> in <module>()
      2 model.initialize()
      3 # run model till termination
----> 4 model.step(till_termination=True)
        global model.step = <bound method AquaCropModel.step of <aquacrop.core.AquaCropModel object at 0x7ff2733f5ba8>>
        global till_termination = undefined

2 frames
/usr/local/lib/python3.6/dist-packages/aquacrop/timestep.py in solution(InitCond=<numba.experimental.jitclass.boxing.InitCondClass object>, ParamStruct=<aquacrop.classes.ParamStructClass object>, ClockStruct=<aquacrop.classes.ClockStructClass object>, weather_step=array([1.7, 6.9, 0.2, 0.0, Timestamp('1980-12-16 00:00:00')], dtype=object), Outputs=<aquacrop.classes.OutputClass object>)
    194 
    195     # 16. Biomass accumulation
--> 196     NewCond = biomass_accumulation(Crop,NewCond,Tr,TrPot_NS,Et0,GrowingSeason)
        NewCond = <numba.experimental.jitclass.boxing.InitCondClass object at 0x7ff26c0e8390>
        global biomass_accumulation = CPUDispatcher(<function biomass_accumulation at 0x7ff27eafb730>)
        Crop = <numba.experimental.jitclass.boxing.CropStruct object at 0x7ff26c0e8670>
        Tr = 0.0
        TrPot_NS = 0.0
        Et0 = 0.0
        GrowingSeason = True
    197 
    198 
itstemo commented 3 years ago

i'm helping @wojeda-mx here, do let us know if know if there was something wrong on our end.

happy to draft a PR if needed, and we're looking forward to keep contributing!

thomasdkelly commented 3 years ago

Hi @wojeda-mx, @itstemo, I have had this issue before and i always fixed it by clipping the the ET0 in the weather file to be at least 0.1mm/day. Not sure whether this has any scientific basis as it was just the minimum ET0 in the Champion Nebraska weather file. I'll ask Tim Foster (the developer of AquaCrop-OS for MATLAB) whether this clipping is sensible and then add it to the prepare_weather function so users dont have to worry about doing it themselves. My guess is that it would be rare for the ET0 in an area to actually be zero so i dont think it is crazy to add a lower limit.

Thanks for spotting the issue and please let me know if you have any other issues or feedback,

thomasdkelly commented 3 years ago

I have added the 0.1 lower limit on ET0 measurements in the latest release 0.1.4.