Seeming issue with orientation.py

[somberostrich]

Version Checks (indicate both or one)

• I have confirmed this bug exists on the lastest release of atlite.

• I have confirmed this bug exists on the current master branch of atlite.

Issue Description

Hi, we've been modelling solar PV output using the latest version of atlite, and we are encountering what appears to be a bug for the 'horizontal' and 'tilted horizontal' tracking options.

For tilted horizontal tracking when the azimuth is set to below 180 degrees (and tilt is set to above 0 degrees), there is a sudden drop in the power output at around noon when the solar zenith approaches 0. This drop is particularly large for the southern hemisphere. In orientation.py, when lines 156 - 160 are commented out:

rotation = np.where(
            logical_and(rotation < 0, azimuth_difference > 0),
            rotation + pi,
            rotation,
        )

the issue appears to disappear and the curves seem to be reasonable.

For horizontal tracking, we also noticed that in orientation.py lines 127 - 129, there might be a bracketing error. It is currently:

surface_azimuth = axis_azimuth + arcsin(
            sin(rotation / sin(surface_slope))
        )

I believe it should be:

surface_azimuth = axis_azimuth + arcsin(
            sin(rotation) / sin(surface_slope)
        )

Reproducible Example

import numpy as np
import pickle
import matplotlib.pyplot as plt
import atlite as at
from datetime import datetime, timedelta
import pytz

class AtlitePV:
    def __init__(self,start_date,end_date,Lon,Lat,tz):
        local_tz = pytz.timezone(tz)
        self.tz=tz
        gmt = pytz.timezone("UTC")

        start_date = local_tz.localize(start_date)
        start_date = start_date.astimezone(gmt)

        end_date = local_tz.localize(end_date)
        end_date = end_date.astimezone(gmt)



        #start_date = start_date.tz_localize(tz).tz_convert(tz='GMT')
        #end_date = end_date.tz_localize(tz).tz_convert(tz='GMT')

        self.cutout = at.Cutout(
            path="ERA-AusSummer-data.nc",
            module="era5",
            x=slice(Lon,Lon+0.75),  # Longitude range
            y=slice(Lat,Lat+0.75),  # Latitude range
            time=slice(str(start_date)[:10],str(end_date)[:10])  # Month or use slice("2025-01-01", "2025-01-31")
        )

        self.cf = None

        self.cutout.prepare(features=["influx", "temperature"])

    def get_solar_power(self,panel:str='CSi',tracking=None,slope:float=0.0,azimuth:float=90.0):
        self.cf = self.cutout.pv(
            panel=panel,
            orientation={"slope": slope, "azimuth": azimuth},
            tracking=tracking,
            shapes=self.cutout.grid,
            #shapes=None,
            per_unit=True
        )

        return self.cf

    def plot_solar_power(self,lon,lat):
        time_index = self.cf.coords['time'].to_index()
        time_local = time_index.tz_localize('UTC').tz_convert('Australia/Melbourne')

        plt.figure(figsize=(12, 4))
        plt.plot(time_local,self.cf)
        plt.title("Solar Capacity Factors")
        plt.ylabel("CF")
        plt.xlabel("Date")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.show()

    def compare_solar_power(self,lon,lat,panel:str='CSi',slope:float=0,azimuth:float=90):
        cf_notrack = self.cutout.pv(
            panel=panel,
            orientation={"slope": slope, "azimuth": azimuth},
            tracking=None,
            shapes=self.cutout.grid,
            # shapes=None,
            per_unit=True
        )

        cf_horizontal = self.cutout.pv(
            panel=panel,
            orientation={"slope": slope, "azimuth": azimuth},
            tracking='horizontal',
            shapes=self.cutout.grid,
            # shapes=None,
            per_unit=True
        )

        cf_horizontal_tilt = self.cutout.pv(
            panel=panel,
            orientation={"slope": slope, "azimuth": azimuth},
            tracking='tilted_horizontal',
            shapes=self.cutout.grid,
            # shapes=None,
            per_unit=True
        )

        cf_vertical = self.cutout.pv(
            panel=panel,
            orientation={"slope": slope, "azimuth": azimuth},
            tracking='vertical',
            shapes=self.cutout.grid,
            # shapes=None,
            per_unit=True
        )

        cf_dual = self.cutout.pv(
            panel=panel,
            orientation={"slope": slope, "azimuth": azimuth},
            tracking='dual',
            shapes=self.cutout.grid,
            # shapes=None,
            per_unit=True
        )

        time_index = self.cf.coords['time'].to_index()
        #time_index = cf_notrack.coords['time'].to_index()
        time_local = time_index.tz_localize('UTC').tz_convert(self.tz)

        fig, ax = plt.subplots(figsize=(12, 4))
        #plt.figure(figsize=(12, 4))
        ax.plot(time_local,np.mean(cf_dual, axis=1),color='black',label="dual tracking")
        ax.plot(time_local, np.mean(cf_vertical, axis=1),color='red',label="vertical tracking")
        ax.plot(time_local, np.mean(cf_horizontal_tilt, axis=1),color='orange',label="tilted horizontal tracking")
        ax.plot(time_local, np.mean(cf_horizontal, axis=1),color='green',label="horizontal tracking")
        ax.plot(time_local, np.mean(cf_notrack, axis=1),color='blue',label="no tracking")

        ax.set_title(f"Atlite output at lat: {lat}, and lon: {lon}")
        ax.set_ylabel("CF")
        ax.set_xlabel("Date")
        ax.grid(True)
        ax.legend()
        plt.show()

        with open(f"track_comp_{Lat}_{Lon}.pkl", "wb") as f:
            pickle.dump(fig, f)

    def plot_tilted(self, lon, lat, panel: str = 'CSi', slope: float = 0, azimuth: float = 90):

        cf_horizontal_tilt = self.cutout.pv(
            panel=panel,
            orientation={"slope": slope, "azimuth": azimuth},
            tracking='tilted_horizontal',
            shapes=self.cutout.grid,
            # shapes=None,
            per_unit=True
        )


        time_index = self.cf.coords['time'].to_index()
        # time_index = cf_notrack.coords['time'].to_index()
        time_local = time_index.tz_localize('UTC').tz_convert(self.tz)

        plt.figure(figsize=(12, 4))
        plt.plot(time_local, np.mean(cf_horizontal_tilt, axis=1), color='orange',
                label="tilted horizontal tracking")

        plt.title(f"Atlite output at lat: {lat}, and lon: {lon}")
        plt.ylabel("CF")
        plt.xlabel("Date")
        plt.grid(True)
        plt.tight_layout()
        plt.show()



if __name__ == "__main__":
    start_date=datetime(2022, 1, 1,hour=1)
    end_date=datetime(2022, 1, 5,hour=23)

    tz = 'Australia/Perth'
    Lon = 115.7
    Lat = -32.2
    azimuth = 0
    slope = 30


    atlite_pv = AtlitePV(start_date=start_date,end_date=end_date,Lon=Lon,Lat=Lat,tz=tz)
    atlite_pv.get_solar_power(azimuth=azimuth,slope=slope)
    atlite_pv.compare_solar_power(lon=Lon, lat=Lat,azimuth=azimuth,slope=slope)

Expected Behavior

No response

Installed Versions

Replace this line.

[euronion]

Thank you for the detailed description and investigation already. Maybe @Parisra who initially authored that part wants to have a look?

Here's the figure from your nice MWE:

[Parisra]

Hi,

Thanks a lot for the clear explanation and the nice example. You are completely correct about the bracketing error for 1-axis horizontal tracking, and I'm including it in a pull request.

As for the drop in generation for tilted horizontal tracking, the issue is actually the incorrect formula for the azimuth difference. The formula currently is:

azimuth_difference = np.where( azimuth_difference > pi, 2 * pi - azimuth_difference, azimuth_difference )

But it should be:

azimuth_difference = np.where( azimuth_difference > pi , azimuth_difference - 2*pi , azimuth_difference )

All the formulas for titled horizontal tracking are based on a short but very useful report from NREL: Rotation Angle for the Optimum Tracking of One-Axis Trackers. As the report states:

the difference γs - γa is evaluated as the angular displacement between two vectors so that it falls within the range of -180° to +180°. For example, if γs =20° and γa = 210°, γs - γa is evaluated as 20° + 360° - 210° = 170° and not 20° - 210° = -190°.

After this correction, the figure from your example would be:

I have now corrected the formula and included it in the pull request as well.

[jwan0038]

Thanks Johannes and Parisra for getting back to me on this, it's much appreciated! I've checked and that seems to work well now.

Cheers,
James