Latitude above 65 °N and problems with plotting results

[Adem6092]

Hi,

I have data from a weather station in northern Sweden (Kiruna, 67.8 °N) and want to evaluate different methods of calculating potential evaporation, which would make PyET a useful tool. I have never used PyET before, and I am not an advanced coder at all. However, following example 1 (as I did) seemed kind of straightforward. But when using latitudes over 65° (converted to radians) as input, I get the following error prompt:

"RuntimeWarning: invalid value encountered in arccos result = getattr(ufunc, method)(*inputs, **kwargs)"

Does the package not support latitudes above a certain point?

I will provide my code below:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyet

data = pd.read_excel(r"met_data_test2.xlsx", index_col=0, parse_dates=True)

meteo = pd.DataFrame({"time":data.index, "tmean":data.T_air_165, "tmax":data.tmax, "tmin":data.tmin, 
                      "rh":data.RH_air_165, "wind":data.Wind_Speed_450, "rs":data.NR*0.0864})

time, tmean, tmax, tmin, rh, wind, rs = [meteo[col] for col in meteo.columns]

elevation = 555 # meters above sea-level
lat = 65.85572 * np.pi/180  # Latitude of the meteorological station, converting from degrees to radians

pet_df = pyet.calculate_all(tmean, wind, rs, elevation, lat, tmax=tmax,tmin=tmin, rh=rh)

fig, axs = plt.subplots(figsize=(13,4), ncols=2)
pet_df.plot(ax=axs[0])
pet_df.cumsum().plot(ax=axs[1], legend=False)

axs[0].set_ylabel("PET [mm/day]", fontsize=12)
axs[1].set_ylabel("Cumulative PET [mm]", fontsize=12)
axs[0].legend(ncol=6, loc=[0,1.])
for i in (0,1):
    axs[i].set_xlabel("Date", fontsize=12)
plt.show()

Thanks for reading my message.

Kind regards,
Adrian

[Adem6092]

I solved it some time ago but i figured that i would share if anyone else would face the same issue - though it is unlikely.

Kiruna (67.8° N) is above the Arctic Circle so the sun doesn't rise for some parts of the winter and doesn't set in some parts of the summer. So i updated the function "sunset_angle" in the source code so that argument to arccos is -1<arg<1:

def sunset_angle(sol_dec, lat):
    if isinstance(lat, DataArray):
        lat = lat.expand_dims(dim={"time": sol_dec.index}, axis=0)
        return arccos(-tan(sol_dec.values) * tan(lat).T).T
    else:
    
        # Calculate the argument of arccos for each element in sol_dec
        arg = -tan(sol_dec) * tan(lat)
        
        # Clamp the argument to the range [-1, 1] for each element
        arg = minimum(maximum(arg, -1), 1)
        
        # Calculate the sunset angle using the clamped argument
        return arccos(arg)

[mvremec]

Dear @Adem6092,

Apologies for the delayed response! We have addressed this issue in the latest version of PyEt 1.3.1. Please update to the most recent version and give it another try. And thank you for your suggestion! We fixed it in a similar way!

If you're still encountering problems after the update, please let me know. I'm here to help!

Cheers,
Matevz

[Adem6092]

Thanks!

The explanation of PM asce on your site probably needs to be updated as it looks like PM now, but according to what I read in the literature it is more similar to PM FAO56 (with less standardized values). The calculation seems right though.

Do you have an idea of which method best estimates potential evaporation from a bare soil? The soil at the site I am working with consists of gravel (or larger grain sizes) and has no vegetation.

[mvremec]

Thank you for identifying the discrepancies in the ASCE documentation. I will update it in the next merge. Regarding your question about estimating potential evaporation from bare soil:

The FAO-56 guidelines recommend using the FAO-56 methods combined with a soil evaporation reduction coefficient. I suggest reviewing the FAO-56 guidelines, specifically pages 142 to 146, for detailed information on how to apply these methods effectively.

[Adem6092]

Thanks, but isn't that en calculation of the acutal evapo(tranpi)ration where they take into account the amount of water in the uppermost "layer"?