Finding the sunniest time for a winter walk
January 6, 2026
My parents live in the French Alps and go on their favorite walk around their village almost daily. During winter, the mountains hide a big part of the sun, as I explained in my previous post. My mom doesn't like the cold and loves the sun. Let’s try to find the optimal time to maximize sun exposure during her walk.
Conveniently, I already know how to import GPS traces from my Garmin watch (see previous posts like this one). Yesterday I went on my parents' usual walk to record the trace:
My parents' usual walk. In many places, the mountains block the sun.
I also know how to compute the horizon using topographic data. All I need to do is put it all together and compute, for a given start time, how much sun there will be on the walk.
Computing the horizon from a given location is relatively expensive. Luckily, the horizon is time independent, it only depends on the location and the azimuth angle (unless a mountain suddenly starts moving, which is unlikely). We can thus precompute the horizon and store it on a grid parameterized by the azimuth angle and the walk duration from the start (the latter is a parameterization of the position along the walk). This is what I obtained:
Precomputed horizon against the azimuth angle (deg) and the position along the walk, parameterized by the walk duration from the start.
Each horizontal slice corresponds to a position in the walk, and each vertical slice corresponds to a viewing direction. There are a few discontinuities at places, probably due to the discretization and finite step I use for the ray casting. These discontinuities should not alter the results too much, I just want an idea of the optimal walk time.
Now that we have the horizon precomputed along the walk for many azimuth angles, we can evaluate this horizon in the direction of the sun at specific dates and times. The sun exposure time is then the amount of time when the sun is above the horizon.
For example, this is what I obtain for today, assuming clear-sky conditions:
Predicted sun exposure during the walk against the starting time.
Anytime between 10am and 12:30pm should provide more than 20 minutes of sun! After that it will feel colder.
Below is what this model predicts for the upcoming weeks:
Predicted sun exposure during the walk in the next few weeks.
Days are getting longer, so is sun exposure! It seems like the walk will be fully in the sun (ignoring trees and buildings and clouds) starting the end of this month, if one starts around noon. On Valentine's day, there will be no need to organize for the perfect sunny walk, since the model suggests that the sun exposure will be maximized whether starting at 10am or 4pm, or anytime between these times. I hope they will go!
The plots have been produced with the following code.
- Plotting the trace:
Show code
#!/usr/bin/env python
import argparse
import json
import matplotlib.pyplot as plt
import numpy as np
from pyproj import Transformer
import rasterio
def main():
parser = argparse.ArgumentParser()
parser.add_argument('data_walk', type=str)
parser.add_argument('dem_path', type=str)
parser.add_argument('--out', type=str, default=None)
args = parser.parse_args()
with open(args.data_walk, 'r') as f:
walk = json.load(f)
lat = np.array(walk['lat'])
lon = np.array(walk['lon'])
with rasterio.open(args.dem_path) as ds:
elevation = ds.read(1)
transform = ds.transform
crs = ds.crs
bounds = ds.bounds
if crs != 'EPSG:32632':
raise ValueError(f'Expected UTM format.')
tf = Transformer.from_crs("EPSG:4326", crs, always_xy=True)
x, y = tf.transform(lon, lat)
fig, ax = plt.subplots(figsize=(6, 3.2))
ax.set_axis_off()
extent = [bounds.left, bounds.right, bounds.bottom, bounds.top]
im = ax.imshow(elevation, cmap='terrain', extent=extent, origin="upper")
ax.plot(x, y, color='r')
ax.set_xlim(292509, 300226)
ax.set_ylim(5098571, 5104322)
cbar = fig.colorbar(im, ax=ax, fraction=0.04, pad=0.02)
cbar.set_label("elevation (m)")
plt.tight_layout()
if args.out is None:
plt.show()
else:
plt.savefig(args.out)
if __name__ == '__main__':
main()
- Precomputing the horizon:
Show code
#!/usr/bin/env python
import argparse
import json
import numpy as np
import pickle
from pyproj import Transformer
import rasterio
from tqdm import tqdm
def find_horizon(ds, elevation, x0, y0, az_deg):
max_d = 50_000 # m
dx, dy = ds.res
step = max(dx, dy) # m
max_steps = int(max_d / step)
theta_max = -np.inf
az = np.deg2rad(az_deg)
ux = np.sin(az)
uy = np.cos(az)
row, col = ds.index(x0, y0)
z0 = elevation[row, col]
# shoot a ray.
for j in range(1, max_steps):
d = j * step
if d >= max_d:
break
x = x0 + ux * d
y = y0 + uy * d
row, col = ds.index(x, y)
if row < 0 or row >= elevation.shape[0]:
break
if col < 0 or col >= elevation.shape[1]:
break
z = elevation[row, col]
theta = np.atan2(z - z0, d)
theta_max = max(theta_max, theta)
horizon_deg = np.rad2deg(theta_max)
return horizon_deg
def main():
parser = argparse.ArgumentParser()
parser.add_argument('data_walk', type=str)
parser.add_argument('dem_path', type=str)
parser.add_argument('--out', type=str, default='horizons.pickle')
args = parser.parse_args()
with open(args.data_walk, 'r') as f:
walk = json.load(f)
lats = []
lons = []
durs = []
# cleanup
for lat, lon, dur in zip(walk['lat'],
walk['lon'],
walk['duration']):
if lat is None or lon is None or dur is None:
continue
lats.append(lat)
lons.append(lon)
durs.append(dur)
lats = np.array(lats)
lons = np.array(lons)
durs = np.array(durs)
# remove duplicates (pauses)
durs_u, idx_u = np.unique(durs, return_index=True)
durs, lats, lons = durs_u, lats[idx_u], lons[idx_u]
# resample every 20s to reduce the number of raycasts
dt = 20
T = durs[-1]
t = np.linspace(0, T, int(T/dt))
lats = np.interp(t, durs, lats)
lons = np.interp(t, durs, lons)
# compute horizons for azimuth from 45 to 315 degrees
# this should cover the whole year in this region
azimuths = np.arange(45, 316, 1)
horizons = np.empty((len(t), len(azimuths)), dtype=float)
with rasterio.open(args.dem_path) as ds:
crs = ds.crs
if crs != 'EPSG:32632':
raise ValueError(f'Expected UTM format.')
tf = Transformer.from_crs("EPSG:4326", ds.crs, always_xy=True)
elevation = ds.read(1)
for iwalk, (lat, lon) in tqdm(enumerate(zip(lats, lons)), total=len(lats)):
x, y = tf.transform(lon, lat)
for iaz, az in enumerate(azimuths):
h = find_horizon(ds, elevation, x, y, az)
horizons[iwalk, iaz] = h
data = {
't': t,
'lat': lats,
'lon': lons,
'az': azimuths,
'horizons': horizons
}
with open(args.out, 'wb') as f:
pickle.dump(data, f, protocol=pickle.HIGHEST_PROTOCOL)
if __name__ == '__main__':
main()
- Plot the horizon on grid:
Show code
#!/usr/bin/env python
import argparse
import matplotlib.pyplot as plt
import pickle
def main():
parser = argparse.ArgumentParser()
parser.add_argument('data', type=str)
parser.add_argument('--out', type=str, default=None)
args = parser.parse_args()
with open(args.data, 'rb') as f:
data = pickle.load(f)
horizons = data['horizons']
lats = data['lat']
lons = data['lon']
azimuths = data['az']
durs = data['t']
fig, ax = plt.subplots()
im = ax.imshow(horizons,
extent=[azimuths[0], azimuths[-1],
durs[0]/60, durs[-1]/60],
origin='lower',
aspect='auto')
cbar = fig.colorbar(im, ax=ax, fraction=0.04, pad=0.02)
cbar.set_label("horizon (deg)")
ax.set_xlabel('azimuth (deg)')
ax.set_ylabel('walk duration (min)')
plt.tight_layout()
if args.out is None:
plt.show()
else:
plt.savefig(args.out)
if __name__ == '__main__':
main()
- Compute the sun exposure along the walk:
Show code
#!/usr/bin/env python
import argparse
import datetime
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pickle
sun_radius = 0.266 # degrees
def day_of_year(d: datetime.date) -> int:
return d.timetuple().tm_yday
def refraction_correction_deg(elevation_deg):
"""
Atmospheric refraction correction in degrees.
Valid for elevations > ~ -1 deg.
"""
h = elevation_deg
if h < -1:
return 0.0
return (1.02 / np.tan(np.deg2rad(h + 10.3 / (h + 5.11)))) / 60.0
def solar_position_noaa(lat_deg, lon_deg, when_local, tz_offset_hours):
"""
Approx solar azimuth/elevation (degrees) using NOAA-style formulas.
lat_deg, lon_deg: observer location (deg), lon positive east.
when_local: naive datetime in LOCAL time (no tzinfo)
tz_offset_hours: e.g. +1 for CET, +2 for CEST
Returns: (azimuth_deg, elevation_deg)
azimuth: degrees clockwise from North (0=N, 90=E)
elevation: degrees above horizon
"""
# --- constants / conversions
deg2rad = np.deg2rad
rad2deg = np.rad2deg
lat = deg2rad(lat_deg)
# Day of year
n = day_of_year(when_local.date())
# Fractional year (gamma) in radians
# hour in local time
hour = when_local.hour + when_local.minute/60 + when_local.second/3600
gamma = 2.0*np.pi/365.0 * (n - 1 + (hour - 12)/24.0)
# Equation of time (minutes)
eqtime = 229.18 * (
0.000075
+ 0.001868*np.cos(gamma)
- 0.032077*np.sin(gamma)
- 0.014615*np.cos(2*gamma)
- 0.040849*np.sin(2*gamma)
)
# Solar declination (radians)
decl = (
0.006918
- 0.399912*np.cos(gamma)
+ 0.070257*np.sin(gamma)
- 0.006758*np.cos(2*gamma)
+ 0.000907*np.sin(2*gamma)
- 0.002697*np.cos(3*gamma)
+ 0.00148*np.sin(3*gamma)
)
# Time offset (minutes)
# lon_deg positive east; tz_offset_hours hours east of UTC
time_offset = eqtime + 4.0*lon_deg - 60.0*tz_offset_hours
# True solar time (minutes)
tst = (hour*60.0 + time_offset) % 1440.0
# Hour angle (degrees)
ha_deg = tst/4.0 - 180.0
ha = deg2rad(ha_deg)
# Solar zenith angle
cos_zenith = np.sin(lat)*np.sin(decl) + np.cos(lat)*np.cos(decl)*np.cos(ha)
cos_zenith = np.clip(cos_zenith, -1.0, 1.0)
zenith = np.arccos(cos_zenith)
elevation = 90.0 - rad2deg(zenith)
# Solar azimuth (NOAA convention)
# Compute azimuth from north, clockwise
sin_az = -np.sin(ha) * np.cos(decl) / np.sin(zenith)
cos_az = (np.sin(decl) - np.sin(lat)*np.cos(zenith)) / (np.cos(lat)*np.sin(zenith))
az = np.arctan2(sin_az, cos_az) # radians, range [-pi, pi]
az_deg = (rad2deg(az) + 360.0) % 360.0
return az_deg, elevation
def find_horizon(ds, elevation, x0, y0, az_deg):
max_d = 20_000 # m
dx, dy = ds.res
step = max(dx, dy) # m
max_steps = int(max_d / step)
theta_max = -np.inf
az = np.deg2rad(az_deg)
ux = np.sin(az)
uy = np.cos(az)
row, col = ds.index(x0, y0)
z0 = elevation[row, col]
# shoot a ray.
for j in range(1, max_steps):
d = j * step
if d >= max_d:
break
x = x0 + ux * d
y = y0 + uy * d
row, col = ds.index(x, y)
if row < 0 or row >= elevation.shape[0]:
break
if col < 0 or col >= elevation.shape[1]:
break
z = elevation[row, col]
theta = np.atan2(z - z0, d)
theta_max = max(theta_max, theta)
horizon_deg = np.rad2deg(theta_max)
return horizon_deg
def compute_solar_time(w_lat, w_lon, w_duration, horizon, tstart, tz_offset_hours=1):
time_sun = 0
time_tot = 0
for lat, lon, dur, dtsec in zip(w_lat[1:],
w_lon[1:],
w_duration[1:],
np.diff(w_duration)):
# get the sun position
t = tstart + datetime.timedelta(seconds=dur)
az, el = solar_position_noaa(lat, lon, t, tz_offset_hours)
el += refraction_correction_deg(el)
# get horizon
hor = horizon(az, dur)
if el + sun_radius >= hor:
time_sun += dtsec
time_tot += dtsec
return time_sun, time_tot
def parse_date(s: str) -> datetime.date:
# accept YYYY-MM-DD
try:
return datetime.date.fromisoformat(s)
except ValueError as e:
raise argparse.ArgumentTypeError("Expected --date in YYYY-MM-DD format") from e
def main():
parser = argparse.ArgumentParser()
parser.add_argument('data', type=str)
parser.add_argument('--date', type=parse_date, default=["2026-01-06"],
nargs='+', help="Local date in YYYY-MM-DD")
parser.add_argument('--start', type=str, default="08:00",
help="Start time (HH:MM) local time")
parser.add_argument('--end', type=str, default="18:00",
help="End time (HH:MM) local time")
parser.add_argument('--tz', type=float, default=1,
help="Timezone offset hours (e.g. 1 for CET, 2 for CEST)")
parser.add_argument('--out', type=str, default=None)
args = parser.parse_args()
with open(args.data, 'rb') as f:
data = pickle.load(f)
horizons = data['horizons']
lats = data['lat']
lons = data['lon']
azimuths = data['az']
durs = data['t']
def horizon(az, t):
az -= azimuths[0]
t -= durs[0]
daz = azimuths[1] - azimuths[0]
ddur = durs[1] - durs[0]
im = min(int(t / ddur), len(durs) - 1)
jm = min(int(az / daz), len(azimuths)-1)
ip = min(im+1, len(durs)-1)
jp = min(jm+1, len(azimuths)-1)
h00 = horizons[im,jm]
h01 = horizons[im,jp]
h10 = horizons[ip,jm]
h11 = horizons[ip,jp]
l0 = (t - im * ddur) / ddur
l1 = (az - jm * daz) / daz
l0 = 1 - l0
l1 = 1 - l1
h0 = h00 * l0 + h10 * (1 - l0)
h1 = h01 * l0 + h11 * (1 - l0)
return h0 * l1 + h1 * (1-l1)
dates = args.date
def parse_hhmm(s: str) -> tuple[int, int]:
try:
hh, mm = s.split(":")
return int(hh), int(mm)
except Exception as e:
raise argparse.ArgumentTypeError("Expected HH:MM") from e
sh, sm = parse_hhmm(args.start)
eh, em = parse_hhmm(args.end)
fig, ax = plt.subplots()
for date in dates:
start_time = []
sun_duration_min = []
tstart = datetime.datetime(date.year, date.month, date.day, sh, sm, 0)
tend = datetime.datetime(date.year, date.month, date.day, eh, em, 0)
steps = int((tend - tstart).total_seconds() / 300) # every 5 min
dt = (tend - tstart) / steps
for i in range(steps + 1):
t = tstart + i * dt
tsun, ttot = compute_solar_time(lats, lons, durs, horizon, t)
start_time.append(t)
sun_duration_min.append(tsun/60)
# minutes since midnight
x = np.array([
t.hour * 60 + t.minute + t.second / 60
for t in start_time
])
ax.plot(x, sun_duration_min, label=tstart.strftime("%B %d, %Y"))
ax.set_xlabel("Time of day")
ax.set_xticks(np.arange(sh*60, eh*60 + 1, 60))
ax.set_xticklabels([f"{h:02d}:00" for h in range(sh, eh+1)])
ax.set_ylabel('sun duration (min)')
ax.legend()
ax.set_xlim(sh*60, eh*60+1)
plt.tight_layout()
if args.out is None:
plt.show()
else:
plt.savefig(args.out)
if __name__ == '__main__':
main()