pvlib using pvgis horizontal data estimates significantly different east and west production

Question

Using pure east and west (i.e. 90 and 270 degrees, respectively) I would assume a similar yearly production for the same location. However, I'm constantly calculating around 100 kWh per kW of panel installed larger production on east facing panels (10%)

latitude, longitude = 55.70018558891628, 12.425528151165556 #somewhre in Copenhagen start_date, end_date = '2015-01-01', '2015-12-31 23:00' surface_tilt = 30

Function I use to get horizontal data in the correct timezone:

def get_pvgis_horizontal_surface(lat, lon, start_date, end_date):
# Define the date range
start = pd.Timestamp(start_date)
end = pd.Timestamp(end_date)
# Retrieve two years of irradiance data from PVGIS:
pvgis_data, meta, _ = pvlib.iotools.get_pvgis_hourly(  
latitude=lat, longitude=lon, start=start_date, end=end_date)  
# Round down to the start of the hour
pvgis_data.index = pvgis_data.index.floor('H')
# extract timezone from index
tz = pvgis_data.index.tz
# convert to local time (tz2)
pvgis_data.index = pvgis_data.index.tz_convert(tz2)
# convert all minutes in index to 00

# Ensure the time series starts from 00:00 of the `start_date`
idx_start = pd.Timestamp(start_date, tz=tz2)
if pvgis_data.index[0] != idx_start:
    # Using loc to get the row corresponding to the first timestamp
    # and then assigning it to the desired timestamp.
    pvgis_data.loc[idx_start] = pvgis_data.loc[pvgis_data.index[0]]
    # Reordering the index to maintain chronological order
    pvgis_data = pvgis_data.sort_index()

# Filter out dates that fall in the next year
end_year = int(end_date.split('-')[0])
pvgis_data = pvgis_data[pvgis_data.index.year <= end_year]

return pvgis_data

ghi.sum() here is 993674.51 (calculated as summed poa_direct and poa_sky_diffuse)

This is the function I use to calculate production from 1 kW of a panel

def calculate_yearly_solar_production(lat, lon, surface_tilt, surface_azimuth, ghi, temp_air, wind_speed):
# create a time range based on the index of one of the inputs
time = ghi.index
site_location = pvlib.location.Location(lat, lon, tz=time.tz)
solpos = site_location.get_solarposition(time)

# calculate irradiance components using the Erbs model
dni_extra = pvlib.irradiance.get_extra_radiation(time, epoch_year=int(time.year[0]))
#
irradianceData = pvlib.irradiance.erbs(ghi, solpos['zenith'], time)
dni = irradianceData['dni'] # does not work well!!!!!
dhi = irradianceData['dhi']
poa_irradiance = pvlib.irradiance.get_total_irradiance(
    surface_tilt, surface_azimuth, solpos['zenith'], solpos['azimuth'],
    dni, ghi, dhi, dni_extra=dni_extra, model='haydavies')
# print(poa_irradiance['poa_global'].describe())

g_poa_effective = poa_irradiance['poa_global'] 
# The temperature coefficient in units of 1/C. Typically -0.002 to -0.005 per degree C.
gamma_pdc = -0.004 # Monocrystalline panels usually have a temperature coefficient between -0.3%/°C and -0.5%/°C.
pdc0=1 # Nameplate DC rating, 1 kW

# define the module parameters
params = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS['pvsyst']['insulated']

# calculate the cell temperature
temperature_cell = pvsyst_cell(poa_irradiance['poa_global'], temp_air, wind_speed, **params)
# print(temperature_cell.describe())

# this is dc power output
dc = pvsystem.pvwatts_dc(g_poa_effective, temperature_cell, pdc0, gamma_pdc, temp_ref=25.0)
    
return dc

So until calculating poa_irradiance, all inputs are the same both for east and for west PV.

However, poa_radiance sum gives significantly different values for poa_direct, while diffuse, sky diffuse and ground diffuse components are rather similar. The difference in direct radiation is what causes the discrepancy in total production.

Sum for the whole 2015: poa_direct west = 400139.545655 poa_direct east = 511274.2

What could be the reason for this significant difference?

Running simulations for east and west for the given inclination directly from pvgis gives rather similar results. So I should also get similar.

Inspecting the hourly results, and comparing my results to pvgis ones, winter production seems to be larger in my case, only for the east facing panels. But what could be the cause here? I checked time zones many times. Am I calculating something wrong with the solar position? Finally, comparing hourly results with pvgis, seems that west facing side is almost exactly the same (the yearly sum is within 2%), so it puzzles me even more that there is a significant discrepancy for the side facing east.