How can I calculate the number of sun hours for a given day and location, considering nearby mountains (with Python)?

Question

My goal is to compute the number of hours for a given day that the sun shines on a given location using Python, under the assumption of a clear sky. The problem arises from the search for real estate. I would like to know how much sun I will actually get on some property, such that I do not have to rely on the statements made by the real estate sales person and can judge a property solely based on the address. I am searching in an area with several nearby mountains, which should be considered in the calcualtion. The approach that I would like to use is the following:

1. For a whole year (1.1.2020 till 31.12.2020), with a given temporal resolution (e.g. minutes), compute the altitude and azimuth angles of the sun for the defined location.
2. Find the angular height of nearby obstacles as seen from the location, yielding value pairs of azimuth and altitude angles. This can be trees, buildings or the already mentioned mountains. Let's assume that trees and buildings are negligible and it's mainly the mountains that take away the sun.
3. For each day, check for each time with the given resolution, whether the altitude angle of the sun is higher than the altitude angle of the obstacle, at the azimuthal location of the sun. For each day of the year, I can then know at which times the sun is visible.

Steps 1 and 3 are easy. For step 1, I can use for example the Python module pysolar. Step 2 is more complicated. If I were to stand on a perfectly flat plane that extends far to the horizon, the obstacle altitude would be 0 for all azimuth angles. If there was a mountain nearby, I would need to know the shape of the mountain, as seen from the location. Unfortunately, I do not even know where to start solving step 2 and I do not know how to name this problem, i.e. I do not know how to google for a solution. In the best case, there would be a Python module that does that the calculation for me, e.g. by connecting to topography data based on OpenStreetMap or other services. If such a module does not exist, I would have to manually program an access to topography data and then do some type of grid search - divide the landscape in a fine grid (possibly spherical coordinates), compute the altitude and azimuth angle of grid points as seen from the location (under consideration of the earth curvature and the elevation at the location) and find the maximum for a given azimuth angle. Are there easier ways to do this? Are there Python modules that do this?

Answer 1

Step 2 isn't that hard if you have a height map and I'm sure there are height maps available.

Cast n rays in an evenly spread 360 degree pattern out from your location and visit all map positions every delta meters from your starting position along that ray. If d is the distance of a point to you, h is its height and h0 is your height on the heightmap keep the maximum value of (h-h0)/d along the path (if something is twice as far away it needs to be twice as high to cast the same length shadow).

You can pretty much ignore earth's curvature - its effect on sunlight occlusion is negligible.