That is amazing. Here is MAthematica code to get this:
sp = SunPosition[
DateRange[DateObject[{2023, 1, 1, 9}, "Hour", "Gregorian", -4.`],
DateObject[{2023, 12, 31, 9}, "Hour", "Gregorian", -4.`], 10]];
Show[
ParametricPlot[
QuantityMagnitude[sp[t]], {t, sp["FirstTime"], sp["LastTime"]}],
Graphics[{Orange, PointSize[Medium],
Point[QuantityMagnitude[sp]["Values"]]}],
Frame -> True
]