Horas de luz por latitud (Daytime by latitude)

daytime_by_latitude.py

import pandas as pd
import matplotlib.pyplot as plt
from datetime import datetime
import math
'''Horas de luz por latitud.
Valores para Mendoza, Argentina; año 2024.
Cambiando los valores de diccionario dlatitudes se puede usar para cualquier lugar de la tierra.
Cambiar también loc1, loc2 con keys a destacar de dlatitudes (horas max/min, equinoccios y solsticios).
Genera dos imágenes png; uno con curvas de horas por latitud y otro de isócronas por latitudes.
Las latitudes en dlatitudes deben ser de punto flotante con 4 decimales.

Daytime by latitude.
Fill dlatitudes dictionary with label and latitude of places of interest.
Also change loc1, loc2 with keys to highlight from dlatitudes (max/min hours, equinox and solstice days).
Generates 2 png images; curves of hours per latitude and isochrones per latitude.
Latitudes in dlatitudes must be floats with 4 decimal values.
'''
def daylength(dayOfYear, lat):
    """Taken from:
    https://gist.github.com/anttilipp/ed3ab35258c7636d87de6499475301ce
    but using math instead of nympy.
    
    Computes the length of the day (the time between sunrise and
    sunset) given the day of the year and latitude of the location.
    Function uses the Brock model for the computations.
    For more information see, for example,
    Forsythe et al., "A model comparison for daylength as a
    function of latitude and day of year", Ecological Modelling,
    1995.

    Parameters
    ----------
    dayOfYear : int
        The day of the year. 1 corresponds to 1st of January
        and 365 to 31st December (on a non-leap year).
    lat : float
        Latitude of the location in degrees. Positive values
        for north and negative for south.

    Returns
    -------
    d : float
        Daylength in hours.
    """
    latInRad = math.radians(lat)
    declinationOfEarth = 23.45*math.sin(math.radians(360.0*(283.0+dayOfYear)/365.0))
    if -math.tan(latInRad) * math.tan(math.radians(declinationOfEarth)) <= -1.0:
        return 24.0
    elif -math.tan(latInRad) * math.tan(math.radians(declinationOfEarth)) >= 1.0:
        return 0.0
    else:
        hourAngle = math.degrees(math.acos(-math.tan(latInRad) * math.tan(math.radians(declinationOfEarth))))
        return 2.0*hourAngle/15.0

def year_day_to_date(year, yday):
    return datetime.strptime(f"{year}-{yday}", "%Y-%j").strftime("%Y-%m-%d")

def date_to_year_day(d):
    return datetime.strptime(d, "%Y-%m-%d").strftime("%j")

def get_xtick_labels(year):
    xticks = {}
    xticks[-50.0] = ""
    xticks[0.0] = f"{year}-01-01"
    for m in range(2, 13):
        for d in [1]:
            xticks[float(date_to_year_day(f"{year}-{m}-{d}"))] = f"{year}-{m:02}-{d:02}"
    return xticks

def format_plot(plt, df, equinoxes, lat_key, lat_value):
    for s in equinoxes: #:
        plt.axvline(x = s, color = 'gray', label = year_day_to_date(year,s),  linestyle = 'dashdot',alpha=0.5, zorder=1)
        plt.text(s-23, 9.4, year_day_to_date(year,s))

    plt.axhline(y = df.min()[lat_value], color = 'gray', label = "Min",  linestyle = 'dotted',alpha=0.5, zorder=1)
    plt.text(385, df.min()[lat_value], f'Mín. {lat_key} ({round(df.min()[lat_value], 2)})')
    plt.axhline(y = df.max()[lat_value], color = 'gray', label = "Max",  linestyle = 'dotted',alpha=0.5, zorder=1)
    plt.text(385, df.max()[lat_value], f'Máx. {lat_key} ({round(df.max()[lat_value], 2)})')

    plt.xlabel(None)
    plt.ylabel('Horas de Luz', fontsize = 14)
    plt.title(f'Horas de Luz [Latitud Sur]', fontsize = 18)


# Create a list of data to 
# be represented in x-axis 
year = 2024
year_days = list(range(1,366))

# lavalle -32.72417230633202, -68.59460899793771
# mendoza -32.88980681767816, -68.822876003072
# tunuyan -33.58149627625932, -69.01873380990953
# san rafael -34.616467633671355, -68.33825669564037
# malargüe -35.473122175850946, -69.58419276211484
# Ranquil Norte -36.65897417899151, -69.83031603573339
# City to latitude dictionary.
dlatitudes = {
    #"Usuhaia": -54.7965,
    #"Trelew": -43.2119,
    "Ranquil Norte": -36.6589,
    "Malargüe": -35.4731,
    "San Rafael": -34.6165,
    "Tunuyán": -33.5815,
    "Mendoza": -32.8898,
    "Capilla del Rosario": -32.1442,
    #"Jujuy": -24.1853,
    #"Quito (Ecuador)": -0.1866
    }
# location of interest
loc1 = "Mendoza"
loc2 = "Ranquil Norte"

latlegends = [ f"{k:20} (Lat: {v})" for k,v in dlatitudes.items()]
latitudes = dlatitudes.values()
# day of year for equinoxes and solstices
equinoxes = []
# build X-axe date labels
year_days_season = [year_day_to_date(year, y) for y in year_days]
xticks = get_xtick_labels(year)

# Hours of daylight per latitude - y axe Hours
lat_hours = {"days" : year_days_season}

# build series of hours of daylight per latitude
for lat in latitudes:
    if lat not in lat_hours.keys():
        lat_hours[lat] = []
    for day in year_days:
        dlen = daylength(day, lat)
        lat_hours[lat].append(round(dlen, 4))
        # equinoxes for latitude of interest
        # if lat == dlatitudes[loc1] and round(dlen,2) >= 11.98 and round(dlen,2) <= 12.01:
        #     # add spring/autumn equinoxes
        #     equinoxes.append(day)

# Create a dataframe
df = pd.DataFrame(lat_hours)
# dataframe for a single location
dfmdz = df[[dlatitudes[loc1]]]
# find spring/autumn equinoxes days (day = night = 12.0 hours)
equinoxes = dfmdz.iloc[(dfmdz[dlatitudes[loc1]] - 12.0).abs().argsort()[:2]].index.values.tolist()
# find max, min hours indexes for a location - summer/winter solstices 
equinoxes.append(dfmdz.idxmax()[dlatitudes[loc1]])
equinoxes.append(dfmdz.idxmin()[dlatitudes[loc1]])
equinoxes.sort()

ax = plt.gca()

#use plot() method on the dataframe 
df.plot( x = 'days', ax = ax, figsize=(20, 12) ) 

handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, latlegends, fontsize=11)
format_plot(plt, dfmdz, equinoxes, loc1, dlatitudes[loc1])
#Add X-axe ticks every first day of month
plt.xticks([ x for x in xticks.keys()], xticks.values(), rotation=30)

plt.savefig('lat-hours.png')
print("Done hours by latitude")
plt.cla()

# Isochrones - latitudes with same hours of daylight - Y axe latitude
# ###################################################################
hours_lat = {"days" : year_days_season}

# Lines of equal hours per latitude (isochrones).
# Transpose Latitud to hours series to Hours to Latitude.
# Massage series for plotting to avoid cluttering the graph.
float_filter = [10.4, 11.4, 11.8, 12.4, 12.8, 13.4]
for i,lt in enumerate(latitudes):
    for d,h in enumerate(df[lt]):
        k = int(h)
        
        # do not show some repeated values that draw as a flat line on top
        if int(k) in [10, 11, 12, 13] and round(h,1) not in float_filter:
            continue
        
        # use float for particular values so those hours are better represented in plot.
        if round(h,1) in float_filter:
            k = round(h,1)
        
        # create an initialize list
        if k not in hours_lat.keys():
            hours_lat[k] = []
            for day in year_days:
                hours_lat[k].insert(day, None)
        
        hours_lat[k][d] = lt


dfhl = pd.DataFrame.from_dict(hours_lat, orient='index')
dfhl = dfhl.transpose()

print("Done rev df")

dfhl.plot( x = 'days', ax = ax,figsize=(15, 5), marker='o',markersize=2)
plt.xlabel(None)
plt.ylabel('Latitud Sur', fontsize = 14)
plt.title("Isócronas por Latitudes")

# horizontal lines for latitudes of interest
for loc in [loc1, loc2]:
    plt.axhline(y = dlatitudes[loc], color = 'gray', label = f"{loc} ({round(dlatitudes[loc],1)})", linestyle = 'dotted', alpha=0.5, zorder=1)
# vertical lines for solstices/equinoxes
for s in equinoxes: #:
    plt.axvline(x = s, color = 'dimgray', label = year_day_to_date(year,s), linestyle = 'dashdot',alpha=0.5, zorder=1)

handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper right', bbox_to_anchor=(0.3, -0.15),ncol=3, fontsize='x-small')
plt.xticks([ x for x in xticks.keys()], xticks.values(), rotation=30, fontsize='x-small')

plt.savefig('lat-hours-iso.png',bbox_extra_artists=(lgd,), bbox_inches='tight')

Los valores de latitude para un lugar se pueden obtener fácilmente de Google Maps con click derecho y copiando el primer valor.

Latitude values can be easily obtained from Google Maps by right clicking a place and copying the first value.

Dos valores "extremos" agregados para comparación, Usuhaia Y Quito (Ecuador).
En Quito las horas diurnas son 12.0 durante todo el año.