surface_azimuth is still ~ 15° off of due South

[Casey McGonigle]

Hi all,

Thanks for your continued support in this google group. It has been invaluable to me as I've gotten exposure to pvlib for the first time.

Anyway, a couple weeks ago I posted in this Google Group asking why my optimal surface azimuth was so far from 180° (actually closer to 0°). Kevin Anderson helped me get that resolved by fixing the timezone -- now my optimal surface_azimuth is at least generally South-facing. However, my code is still returning an optimal surface_azimuth that is ~ 15° East of South (for tilt = 25°).

I've tried using multiple cities throughout the US (Berkeley, CA; Rochester, NY; Kansas City, MO; Austin, TX) and they all return an optimal surface_azimuth that is at least 14° East of due South.

Moreover, if I increase the surface_tilt angle (say to vertical, 90°), the optimal surface_azimuth diverges even further from 180°. In Berkeley, that optimal surface_azimuth is 138°; Austin is the worst. It goes all the way to 108°.

Thanks for reading. I really appreciate any help you can offer.

I've copied my code below and attached a couple of visualizations to help with troubleshooting

```

import os

import pandas as pd

import numpy as np

import os

import os.path

import matplotlib.pyplot as plt

import pvlib

from geopy.exc import GeocoderTimedOut 

from geopy.geocoders import Nominatim 

from IPython.display import Image

from timezonefinder import TimezoneFinder

## GET THE LATITUDE & LONGITUDE OF A GIVEN CITY

geolocator = Nominatim(user_agent="____'s app") 

geo = geolocator.geocode("Austin") 

## CHECK THAT CITY IS CORRECT (by Country, State, etc.)

print(geo.address)

# CHECK THE LAT, LON order

print(geo.latitude, geo.longitude)

## CREATE OUR TIMEZONE

tf = TimezoneFinder()

#latitude, longitude = 52.5061, 13.358

tzone = tf.timezone_at(lng=geo.longitude, lat=geo.latitude)

print(tzone)

## SELECT THE YEAR & TIMES YOU'D LIKE TO MODEL OFF

YEAR = 2019

STARTDATE = '%d-01-01T00:00:00' % YEAR

ENDDATE = '%d-12-31T23:59:59' % YEAR

TIMES = pd.date_range(start=STARTDATE, end=ENDDATE, freq='H', tz = tzone)

## ACCESS THE NREL API TO EXTRACT WEATHER DATA

NREL_API_KEY = os.getenv('NREL_API_KEY', 'DEMO_KEY')

## fill in the blank with your email address

EMAIL = os.getenv('EMAIL', 'cmc...@berkeley.edu')

##NEED TO COMMENT OUT THIS LINE BELOW -- if you call it too many times within an hour, it will break your code

header, data = pvlib.iotools.get_psm3(geo.latitude, geo.longitude, NREL_API_KEY, EMAIL)

## SELECT THE PVLIB PANEL & INTERVTER YOU'D LIKE TO USE

## CAN ALSO CHOOSE FROM SANDIA LABS' DATASET OF PANELS & INVERTERS (check out the function)

## WE CHOSE THE CECMods because they highlighted the ones that were BIPV

INVERTERS = pvlib.pvsystem.retrieve_sam('CECInverter')

INVERTER_10K = INVERTERS['SMA_America__SB10000TL_US__240V_']

CECMODS = pvlib.pvsystem.retrieve_sam('CECMod')

## select the panel you'd like to use (note: this is a BIPV panel)

CECMOD_MONO = CECMODS['Pevafersa_America_IP_235_GG']

## CREATING AN ARRAY TO ITERATE THROUGH IN ORDER TO TEST DIFFERENT SURFACE_AZIMUTHS

heading_array = np.arange(0, 361, 2)

heading_DF = pd.DataFrame(heading_array).rename(columns = {0: "Heading"})

heading_DF.head()

# geo IS AN OBJECT (the given city) CREATED ABOVE

LATITUDE, LONGITUDE = geo.latitude, geo.longitude

# data IS AN OBJECT (the weather patterns) CREATED ABOVE

# TIMES IS ALSO CREATED ABOVE, AND REPRESENTS TIME

data.index = TIMES

dni = data.DNI.values

ghi = data.GHI.values

dhi = data.DHI.values

surface_albedo = data['Surface Albedo'].values

temp_air = data.Temperature.values

dni_extra = pvlib.irradiance.get_extra_radiation(TIMES).values

# GET SOLAR POSITION

sp = pvlib.solarposition.get_solarposition(TIMES, LATITUDE, LONGITUDE)

solar_zenith = sp.apparent_zenith.values

solar_azimuth = sp.azimuth.values

# CREATING THE ARRY TO STORE THE DAILY ENERGY OUTPUT BY SOLAR AZIMUTH

e_by_az = []

# IDEAL ANGLE IN NORTHERN HEMISPHERE IS ~25; vertical is ~90

surface_tilt = 25

for heading in heading_DF["Heading"]:

    surface_azimuth = heading

    poa_sky_diffuse = pvlib.irradiance.get_sky_diffuse(

        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,

        dni, ghi, dhi, dni_extra=dni_extra, model='haydavies')

    # calculate the angle of incidence using the surface_azimuth and (hardcoded) surface_tilt

    aoi = pvlib.irradiance.aoi(

        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth)

    # https://pvlib-python.readthedocs.io/en/stable/generated/pvlib.irradiance.aoi.html

    # https://pvlib-python.readthedocs.io/en/stable/generated/pvlib.pvsystem.PVSystem.html

    poa_ground_diffuse = pvlib.irradiance.get_ground_diffuse(

        surface_tilt, ghi, albedo=surface_albedo)

    poa = pvlib.irradiance.poa_components(

        aoi, dni, poa_sky_diffuse, poa_ground_diffuse)

    poa_direct = poa['poa_direct']

    poa_diffuse = poa['poa_diffuse']

    poa_global = poa['poa_global']

    iam = pvlib.iam.ashrae(aoi)

    effective_irradiance = poa_direct*iam + poa_diffuse

    temp_cell = pvlib.temperature.pvsyst_cell(poa_global, temp_air)

    # THIS IS THE MAGIC

    cecparams = pvlib.pvsystem.calcparams_cec(

        effective_irradiance, temp_cell,

        CECMOD_MONO.alpha_sc, CECMOD_MONO.a_ref,

        CECMOD_MONO.I_L_ref, CECMOD_MONO.I_o_ref,

        CECMOD_MONO.R_sh_ref, CECMOD_MONO.R_s, CECMOD_MONO.Adjust)

    # mpp is the list of energy output by hour for the whole year using a single panel

    mpp = pvlib.pvsystem.max_power_point(*cecparams, method='newton')

    mpp = pd.DataFrame(mpp, index=TIMES)

    first48 = mpp[:48]

    Edaily = mpp.p_mp.resample('D').sum()

    # Edaily is the list of energy output by day for the whole year using a single panel

    Eyearly = sum(Edaily)

    

    e_by_az.append(Eyearly)

    

## LINKING THE Heading (ie. surface_azimuth) AND THE Eyearly (ie. yearly energy output) IN A DF

heading_DF["Eyearly"] = e_by_az

heading_DF.head()

## VISUALIZE ENERGY OUTPUT BY SURFACE_AZIMUTH

## PLEASE SEE THE ATTACHED VISUALIZATIONS

plt.plot(heading_DF["Heading"], heading_DF["Eyearly"])

plt.xlabel("Surface_Azimuth Angle")

plt.ylabel("Yearly Energy Output with tilt @ " + str(surface_tilt))

plt.title("Yearly Energy Output by Solar_Azimuth Angle using surface_tilt = " + str(surface_tilt) + " in Austin, TX");

```

Please see the attached visualizations to understand my output/issue

I really appreciate your help,

Casey

[Kevin Anderson]

Hi Casey,

I think there may be a timezone issue here :)  The TimezoneFinder method you're using returns timezones like "America/Chicago", which include the effect of DST.  However, the data returned by iotools.get_psm3 has timezones like "Etc/GMT+6", which have no DST.  So overwriting the non-DST index to one with DST will introduce a 1-hour error for the part of the year when DST is active.

But stepping back from that -- why bother overwriting the data index at all?  Is it to collapse the various years comprising the TMY into one continuous year of data?  If so, you might consider just editing the years (leaving the timezone info intact) instead of creating a whole new TIMES index:

    data.index = map(lambda dt: dt.replace(year=YEAR), data.index)

After using the above line and removing all the TIMES references, I get these "optimal azimuths":

- Berkeley: 186
- Rochester: 180
- Kansas City: 180
- Austin: 190
- Denver: 174
- Raleigh: 178

Austin being 10 degrees "off" is interesting.  I didn't look into it much, but here's a heatmap of GHI with solar noon overlaid in black; my qualitative impression is that there tends to be more irradiance in the afternoon, which is consistent with optimizing to a slightly western azimuth.  Not sure if this is indicative of some other problem or just a quirk of Austin's climate.

Cheers,
Kevin

--
You received this message because you are subscribed to the Google Groups "pvlib-python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pvlib-python...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pvlib-python/63466955-f173-4ed5-996a-b57104d7ef4fn%40googlegroups.com.

[Casey McGonigle]

Kevin!

Thank you so much, you've been so helpful.

I implemented your "data.index = map(lambda dt: dt.replace(year=YEAR), data.index)" and it's working great for me -- I was able to corroborate your results with surface_tilt = 25°. I never would have found that bug without your help.

In my first email, I specifically mentioned "Moreover, if I increase the surface_tilt angle (say to vertical, 90°), the optimal surface_azimuth diverges even further from 180°" because I'm working to model BIPV windows on the sides of houses (ie. surface_tilt = 90°). That further divergence at surface_tilt = 90° is still apparent after I implemented your solution. Here are the "optimal azimuths" with surface_tilt = 90°:

- Berkeley: 212 (peaks away from 180)
- Rochester: 202 (but flat peak includes 180)
- Kansas City: 172 (flat peak includes 180)
- Austin: 224 (peak away from 180)
- Denver: 162 (flat peak includes 180)
- Raleigh: 166 (flat peak includes 180)

I've also attached graphs of their Wattage by surface_azimuth with surface_tilt = 90°

Does it make sense to you that vertical panels would have significantly different (and further from 180°) optimal azimuths than panels with a more optimal tilt?

Thanks for all your help,

Casey

[Kevin Anderson]

Having never considered vertical panels before, I can only offer speculation, but here it is.  Consider, if instead of tilt=25, you set 90 < tilt < 180, so the panels were tilted over with the front facing somewhat down instead of somewhat up?  If the azimuth is 180, you'll hardly collect any direct light because when the sun is low (morning and evening) it's way off to the side, and in midday it's too high in the sky for the sun to see the front face.  But if the panel were instead facing east or west instead of south (still with tilt>90), it could at least collect some direct irradiance in the morning or afternoon before the sun gets too high or too far south in the sky.  That would produce a double-peaked optimization profile, presumably roughly symmetric around 180.  And finally, if there is a double-peaked profile for tilt > 90 and single-peaked profile for tilt < 90, there must be a broad flat peak for tilts around 90 as the boundary between single- and double-peaked. 

I added a nested loop to your code that iterates across surface_tilts and here is the result showing this more concretely.  Color indicates annual production and dark lines indicate contours of equal production.  Somewhere around tilt=90, the single peak at 180 bifurcates and heads off towards the east and west.  Note that modeling tilt=170 is kind of a misuse of pvlib's models (e.g. this ignores that module temperature would switch to be driven more by rear-side irradiance than front-side) but I wouldn't expect those conceptual errors to change the shape of these curves much.

Perhaps someone else will chime in to comment on whether this intuitive explanation is correct.

Cheers,

Kevin

To view this discussion on the web visit https://groups.google.com/d/msgid/pvlib-python/64ce037b-ac9e-476f-835a-58e0babf087an%40googlegroups.com.

[Casey McGonigle]

Kevin,

Thank you very much, once again! Your explanation (and awesome graph) make intuitive sense to me and your commitment to rapid, in-depth replies makes a big difference to people like me.

Thank you for all your help,

Casey

[Casey McGonigle]

ps is there any chance you can share the code for either or both of these graphs? I'd love to play around with them, but know it would take me a couple hours to get them working.

Thanks again,
Casey

[kevina...@gmail.com]

Thanks for asking good, specific questions! 

I didn't save that code (d'oh), but I've attached a modified version of your script that is more or less what I was using.  The two important bits for plotting are df.pivot_table to fold a series into a matrix and plt.pcolormesh to display it as a heatmap.

Kevin

[Casey McGonigle]

Awesome! I really appreciate all your help,

Casey