TMY to Power Tutorial

This tutorial will walk through the process of going from TMY data to AC power using the SAPM.

Table of contents:

  1. Setup
  2. Load TMY data
  3. Calculate modeling intermediates
  4. DC power using SAPM
  5. AC power using SAPM

This tutorial has been tested against the following package versions:

  • pvlib 0.2.0
  • Python 2.7.10
  • IPython 3.2
  • pandas 0.16.2

It should work with other Python and Pandas versions. It requires pvlib >= 0.2.0 and IPython >= 3.0.

Authors:

  • Will Holmgren (@wholmgren), University of Arizona, July 2015
  • Rob Andrews (@Calama-Consulting), Heliolytics, June 2014

Setup

These are just your standard interactive scientific python imports that you'll get very used to using.


In [1]:
# built-in python modules
import os
import inspect

# scientific python add-ons
import numpy as np
import pandas as pd

# plotting stuff
# first line makes the plots appear in the notebook
%matplotlib inline 
import matplotlib.pyplot as plt
import matplotlib as mpl
# seaborn makes your plots look better
try:
    import seaborn as sns
    sns.set(rc={"figure.figsize": (12, 6)})
    sns.set_color_codes()
except ImportError:
    print('We suggest you install seaborn using conda or pip and rerun this cell')

# finally, we import the pvlib library
import pvlib

Load TMY data

pvlib comes with a couple of TMY files, and we'll use one of them for simplicity. You could also load a file from disk, or specify a url. See this NREL website for a list of TMY files:

http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2005/tmy3/by_state_and_city.html


In [2]:
# Find the absolute file path to your pvlib installation
pvlib_abspath = os.path.dirname(os.path.abspath(inspect.getfile(pvlib)))

# absolute path to a data file
datapath = os.path.join(pvlib_abspath, 'data', '703165TY.csv')

# read tmy data with year values coerced to a single year
tmy_data, meta = pvlib.tmy.readtmy3(datapath, coerce_year=2015)
tmy_data.index.name = 'Time'

# TMY data seems to be given as hourly data with time stamp at the end
# shift the index 30 Minutes back for calculation of sun positions
tmy_data = tmy_data.shift(freq='-30Min')

The file handling above looks complicated because we're trying to account for the many different ways that people will run this notebook on their systems. You can just put a simple string path into the readtmy3 function if you know where the file is.

Let's look at the imported version of the TMY file.


In [3]:
tmy_data.head()


Out[3]:
ETR ETRN GHI GHISource GHIUncertainty DNI DNISource DNIUncertainty DHI DHISource ... AOD AODSource AODUncertainty Alb AlbSource AlbUncertainty Lprecipdepth Lprecipquantity LprecipSource LprecipUncertainty
Time
2015-01-01 00:30:00-09:00 0 0 0 1 0 0 1 0 0 1 ... 0.051 F 8 0.24 F 8 -9900 -9900 ? 0
2015-01-01 01:30:00-09:00 0 0 0 1 0 0 1 0 0 1 ... 0.051 F 8 0.24 F 8 -9900 -9900 ? 0
2015-01-01 02:30:00-09:00 0 0 0 1 0 0 1 0 0 1 ... 0.051 F 8 0.24 F 8 -9900 -9900 ? 0
2015-01-01 03:30:00-09:00 0 0 0 1 0 0 1 0 0 1 ... 0.051 F 8 0.24 F 8 -9900 -9900 ? 0
2015-01-01 04:30:00-09:00 0 0 0 1 0 0 1 0 0 1 ... 0.051 F 8 0.24 F 8 -9900 -9900 ? 0

5 rows × 66 columns

This is a pandas DataFrame object. It has a lot of great properties that are beyond the scope of our tutorials.

Plot the GHI data from the TMY file


In [4]:
tmy_data['GHI'].plot()
plt.ylabel('Irradiance (W/m**2)')


Out[4]:
<matplotlib.text.Text at 0x87a1320>

Calculate modeling intermediates

Before we can calculate power for all times in the TMY file, we will need to calculate:

  • solar position
  • extra terrestrial radiation
  • airmass
  • angle of incidence
  • POA sky and ground diffuse radiation
  • cell and module temperatures

First, define some PV system parameters.


In [5]:
surface_tilt = 30
surface_azimuth = 180 # pvlib uses 0=North, 90=East, 180=South, 270=West convention
albedo = 0.2

# create pvlib Location object based on meta data
sand_point = pvlib.location.Location(meta['latitude'], meta['longitude'], tz='US/Alaska', 
                                     altitude=meta['altitude'], name=meta['Name'].replace('"',''))
print(sand_point)


SAND POINT: latitude=55.317, longitude=-160.517, tz=US/Alaska, altitude=7.0

Solar position

Calculate the solar position for all times in the TMY file.

The default solar position algorithm is based on Reda and Andreas (2004). Our implementation is pretty fast, but you can make it even faster if you install numba and use add method='nrel_numba' to the function call below.


In [6]:
print(tmy_data.index)
solpos = pvlib.solarposition.get_solarposition(tmy_data.index, sand_point)

solpos.plot()


DatetimeIndex(['2015-01-01 00:30:00-09:00', '2015-01-01 01:30:00-09:00',
               '2015-01-01 02:30:00-09:00', '2015-01-01 03:30:00-09:00',
               '2015-01-01 04:30:00-09:00', '2015-01-01 05:30:00-09:00',
               '2015-01-01 06:30:00-09:00', '2015-01-01 07:30:00-09:00',
               '2015-01-01 08:30:00-09:00', '2015-01-01 09:30:00-09:00', 
               ...
               '2015-12-31 14:30:00-09:00', '2015-12-31 15:30:00-09:00',
               '2015-12-31 16:30:00-09:00', '2015-12-31 17:30:00-09:00',
               '2015-12-31 18:30:00-09:00', '2015-12-31 19:30:00-09:00',
               '2015-12-31 20:30:00-09:00', '2015-12-31 21:30:00-09:00',
               '2015-12-31 22:30:00-09:00', '2014-12-31 23:30:00-09:00'],
              dtype='datetime64[ns]', name='Time', length=8760, freq=None, tz='pytz.FixedOffset(-540)')
Out[6]:
<matplotlib.axes._subplots.AxesSubplot at 0x89848d0>

The funny looking jump in the azimuth is just due to the coarse time sampling in the TMY file.

DNI ET

Calculate extra terrestrial radiation. This is needed for many plane of array diffuse irradiance models.


In [7]:
# the extraradiation function returns a simple numpy array
# instead of a nice pandas series. We will change this
# in a future version
dni_extra = pvlib.irradiance.extraradiation(tmy_data.index)
dni_extra = pd.Series(dni_extra, index=tmy_data.index)

dni_extra.plot()
plt.ylabel('Extra terrestrial radiation (W/m**2)')


Out[7]:
<matplotlib.text.Text at 0x8ad1780>

Airmass

Calculate airmass. Lots of model options here, see the atmosphere module tutorial for more details.


In [8]:
airmass = pvlib.atmosphere.relativeairmass(solpos['apparent_zenith'])

airmass.plot()
plt.ylabel('Airmass')


Out[8]:
<matplotlib.text.Text at 0x8a0dd30>

The funny appearance is due to aliasing and setting invalid numbers equal to NaN. Replot just a day or two and you'll see that the numbers are right.

POA sky diffuse

Use the Hay Davies model to calculate the plane of array diffuse sky radiation. See the irradiance module tutorial for comparisons of different models.


In [9]:
tmy_data['DNI']


Out[9]:
Time
2015-01-01 00:30:00-09:00      0
2015-01-01 01:30:00-09:00      0
2015-01-01 02:30:00-09:00      0
2015-01-01 03:30:00-09:00      0
2015-01-01 04:30:00-09:00      0
2015-01-01 05:30:00-09:00      0
2015-01-01 06:30:00-09:00      0
2015-01-01 07:30:00-09:00      0
2015-01-01 08:30:00-09:00      0
2015-01-01 09:30:00-09:00      0
2015-01-01 10:30:00-09:00      0
2015-01-01 11:30:00-09:00      0
2015-01-01 12:30:00-09:00      0
2015-01-01 13:30:00-09:00      0
2015-01-01 14:30:00-09:00      0
2015-01-01 15:30:00-09:00      0
2015-01-01 16:30:00-09:00      0
2015-01-01 17:30:00-09:00      0
2015-01-01 18:30:00-09:00      0
2015-01-01 19:30:00-09:00      0
2015-01-01 20:30:00-09:00      0
2015-01-01 21:30:00-09:00      0
2015-01-01 22:30:00-09:00      0
2015-01-01 23:30:00-09:00      0
2015-01-02 00:30:00-09:00      0
2015-01-02 01:30:00-09:00      0
2015-01-02 02:30:00-09:00      0
2015-01-02 03:30:00-09:00      0
2015-01-02 04:30:00-09:00      0
2015-01-02 05:30:00-09:00      0
                            ... 
2015-12-30 18:30:00-09:00      0
2015-12-30 19:30:00-09:00      0
2015-12-30 20:30:00-09:00      0
2015-12-30 21:30:00-09:00      0
2015-12-30 22:30:00-09:00      0
2015-12-30 23:30:00-09:00      0
2015-12-31 00:30:00-09:00      0
2015-12-31 01:30:00-09:00      0
2015-12-31 02:30:00-09:00      0
2015-12-31 03:30:00-09:00      0
2015-12-31 04:30:00-09:00      0
2015-12-31 05:30:00-09:00      0
2015-12-31 06:30:00-09:00      0
2015-12-31 07:30:00-09:00      0
2015-12-31 08:30:00-09:00      0
2015-12-31 09:30:00-09:00      0
2015-12-31 10:30:00-09:00     57
2015-12-31 11:30:00-09:00    132
2015-12-31 12:30:00-09:00    531
2015-12-31 13:30:00-09:00    583
2015-12-31 14:30:00-09:00    566
2015-12-31 15:30:00-09:00    474
2015-12-31 16:30:00-09:00    267
2015-12-31 17:30:00-09:00     33
2015-12-31 18:30:00-09:00      0
2015-12-31 19:30:00-09:00      0
2015-12-31 20:30:00-09:00      0
2015-12-31 21:30:00-09:00      0
2015-12-31 22:30:00-09:00      0
2014-12-31 23:30:00-09:00      0
Name: DNI, dtype: int64

In [10]:
poa_sky_diffuse = pvlib.irradiance.haydavies(surface_tilt, surface_azimuth,
                                             tmy_data['DHI'], tmy_data['DNI'], dni_extra,
                                             solpos['apparent_zenith'], solpos['azimuth'])

poa_sky_diffuse.plot()
plt.ylabel('Irradiance (W/m**2)')


Out[10]:
<matplotlib.text.Text at 0x918d080>

POA ground diffuse

Calculate ground diffuse. We specified the albedo above. You could have also provided a string to the surface_type keyword argument.


In [11]:
poa_ground_diffuse = pvlib.irradiance.grounddiffuse(surface_tilt, tmy_data['GHI'], albedo=albedo)

poa_ground_diffuse.plot()
plt.ylabel('Irradiance (W/m**2)')


Out[11]:
<matplotlib.text.Text at 0x9101ac8>

AOI

Calculate AOI


In [12]:
aoi = pvlib.irradiance.aoi(surface_tilt, surface_azimuth, solpos['apparent_zenith'], solpos['azimuth'])

aoi.plot()
plt.ylabel('Angle of incidence (deg)')


Out[12]:
<matplotlib.text.Text at 0xa5e5d68>

Note that AOI has values greater than 90 deg. This is ok.

POA total

Calculate POA irradiance


In [13]:
poa_irrad = pvlib.irradiance.globalinplane(aoi, tmy_data['DNI'], poa_sky_diffuse, poa_ground_diffuse)

poa_irrad.plot()
plt.ylabel('Irradiance (W/m**2)')
plt.title('POA Irradiance')


Out[13]:
<matplotlib.text.Text at 0xbfaef98>

Cell and module temperature

Calculate pv cell and module temperature


In [14]:
pvtemps = pvlib.pvsystem.sapm_celltemp(poa_irrad['poa_global'], tmy_data['Wspd'], tmy_data['DryBulb'])

pvtemps.plot()
plt.ylabel('Temperature (C)')


Out[14]:
<matplotlib.text.Text at 0xa6ad128>

DC power using SAPM

Get module data from the web.


In [15]:
sandia_modules = pvlib.pvsystem.retrieve_sam(name='SandiaMod')

Choose a particular module


In [16]:
sandia_module = sandia_modules.Canadian_Solar_CS5P_220M___2009_
sandia_module


Out[16]:
Vintage                                                          2009
Area                                                            1.701
Material                                                         c-Si
Cells_in_Series                                                    96
Parallel_Strings                                                    1
Isco                                                          5.09115
Voco                                                          59.2608
Impo                                                          4.54629
Vmpo                                                          48.3156
Aisc                                                         0.000397
Aimp                                                         0.000181
C0                                                            1.01284
C1                                                         -0.0128398
Bvoco                                                        -0.21696
Mbvoc                                                               0
Bvmpo                                                       -0.235488
Mbvmp                                                               0
N                                                              1.4032
C2                                                           0.279317
C3                                                           -7.24463
A0                                                           0.928385
A1                                                           0.068093
A2                                                         -0.0157738
A3                                                          0.0016606
A4                                                          -6.93e-05
B0                                                                  1
B1                                                          -0.002438
B2                                                          0.0003103
B3                                                         -1.246e-05
B4                                                           2.11e-07
B5                                                          -1.36e-09
DTC                                                                 3
FD                                                                  1
A                                                            -3.40641
B                                                          -0.0842075
C4                                                           0.996446
C5                                                           0.003554
IXO                                                           4.97599
IXXO                                                          3.18803
C6                                                            1.15535
C7                                                          -0.155353
Notes               Source: Sandia National Laboratories Updated 9...
Name: Canadian_Solar_CS5P_220M___2009_, dtype: object

Run the SAPM using the parameters we calculated above.


In [17]:
sapm_out = pvlib.pvsystem.sapm(sandia_module, poa_irrad.poa_direct, poa_irrad.poa_diffuse,
                              pvtemps['temp_cell'], airmass, aoi)
print(sapm_out.head())

sapm_out[['p_mp']].plot()
plt.ylabel('DC Power (W)')


                           i_sc  i_mp  v_oc  v_mp  p_mp  i_x  i_xx  \
Time                                                                 
2015-01-01 09:30:00+00:00     0     0     0     0     0    0     0   
2015-01-01 10:30:00+00:00     0     0     0     0     0    0     0   
2015-01-01 11:30:00+00:00     0     0     0     0     0    0     0   
2015-01-01 12:30:00+00:00     0     0     0     0     0    0     0   
2015-01-01 13:30:00+00:00     0     0     0     0     0    0     0   

                           effective_irradiance  
Time                                             
2015-01-01 09:30:00+00:00                     0  
2015-01-01 10:30:00+00:00                     0  
2015-01-01 11:30:00+00:00                     0  
2015-01-01 12:30:00+00:00                     0  
2015-01-01 13:30:00+00:00                     0  
Out[17]:
<matplotlib.text.Text at 0xc01a588>

DC power using single diode


In [18]:
cec_modules = pvlib.pvsystem.retrieve_sam(name='CECMod')
cec_module = cec_modules.Canadian_Solar_CS5P_220M

In [19]:
photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth = (
    pvlib.pvsystem.calcparams_desoto(poa_irrad.poa_global,
                                 temp_cell=pvtemps['temp_cell'],
                                 alpha_isc=cec_module['alpha_sc'],
                                 module_parameters=cec_module,
                                 EgRef=1.121,
                                 dEgdT=-0.0002677) )

In [20]:
single_diode_out = pvlib.pvsystem.singlediode(cec_module, photocurrent, saturation_current,
                                              resistance_series, resistance_shunt, nNsVth)

In [21]:
single_diode_out[['p_mp']].plot()
plt.ylabel('DC Power (W)')


Out[21]:
<matplotlib.text.Text at 0x106594e0>

AC power using SAPM

Get the inverter database from the web


In [23]:
sapm_inverters = pvlib.pvsystem.retrieve_sam('sandiainverter')

Choose a particular inverter


In [24]:
sapm_inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_']
sapm_inverter


Out[24]:
Vac          208.000000
Paco         250.000000
Pdco         259.522050
Vdco          40.242603
Pso            1.771614
C0            -0.000025
C1            -0.000090
C2             0.000669
C3            -0.018900
Pnt            0.020000
Vdcmax        65.000000
Idcmax        10.000000
Mppt_low      20.000000
Mppt_high     50.000000
Name: ABB__MICRO_0_25_I_OUTD_US_208_208V__CEC_2014_, dtype: float64

In [25]:
p_acs = pd.DataFrame()
p_acs['sapm'] = pvlib.pvsystem.snlinverter(sapm_inverter, sapm_out.v_mp, sapm_out.p_mp)
p_acs['sd'] = pvlib.pvsystem.snlinverter(sapm_inverter, single_diode_out.v_mp, single_diode_out.p_mp)

p_acs.plot()
plt.ylabel('AC Power (W)')


Out[25]:
<matplotlib.text.Text at 0x1075c7f0>

In [26]:
diff = p_acs['sapm'] - p_acs['sd']
diff.plot()
plt.ylabel('SAPM - SD Power (W)')


Out[26]:
<matplotlib.text.Text at 0x1095db70>

Plot just a few days.


In [27]:
p_acs['2015-07-05':'2015-07-06'].plot()


Out[27]:
<matplotlib.axes._subplots.AxesSubplot at 0x1098def0>

Some statistics on the AC power


In [28]:
p_acs.describe()


Out[28]:
sapm sd
count 8760.000000 4620.000000
mean 23.092072 45.975969
std 42.338813 50.010571
min -0.020000 -0.020000
25% -0.020000 8.625879
50% -0.020000 28.616357
75% 27.946162 63.817331
max 211.292680 213.052674

In [29]:
p_acs.sum()


Out[29]:
sapm    202286.553521
sd      212408.976960
dtype: float64

In [30]:
# create data for a y=x line
p_ac_max = p_acs.max().max()
yxline = np.arange(0, p_ac_max)

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111, aspect='equal')
sc = ax.scatter(p_acs['sd'], p_acs['sapm'], c=poa_irrad.poa_global, alpha=1, cmap=mpl.cm.YlGnBu_r)  
ax.plot(yxline, yxline, 'r', linewidth=3)
ax.set_xlim(0, None)
ax.set_ylim(0, None)
ax.set_xlabel('Single Diode model')
ax.set_ylabel('Sandia model')
fig.colorbar(sc, label='POA Global (W/m**2)')


Out[30]:
<matplotlib.colorbar.Colorbar at 0x10b12d30>

We can change the value of color value c to see the sensitivity of model accuracy to measured meterological conditions. It can be useful to define a simple plotting function for this kind of exploratory analysis.


In [31]:
def sapm_sd_scatter(c_data, label=None, **kwargs):
    """Display a scatter plot of SAPM p_ac vs. single diode p_ac.
    
    You need to re-execute this cell if you re-run the p_ac calculation.
    
    Parameters
    ----------
    c_data : array-like
        Determines the color of each point on the scatter plot.
        Must be same length as p_acs.
        
    kwargs passed to ``scatter``.
    
    Returns
    -------
    tuple of fig, ax objects
    """
    
    fig = plt.figure(figsize=(12,12))
    ax = fig.add_subplot(111, aspect='equal')
    sc = ax.scatter(p_acs['sd'], p_acs['sapm'], c=c_data, alpha=1, cmap=mpl.cm.YlGnBu_r, **kwargs)  
    ax.plot(yxline, yxline, 'r', linewidth=3)
    ax.set_xlim(0, None)
    ax.set_ylim(0, None)
    ax.set_xlabel('Single diode model power (W)')
    ax.set_ylabel('Sandia model power (W)')
    fig.colorbar(sc, label='{}'.format(label), shrink=0.75)
    
    return fig, ax

In [32]:
sapm_sd_scatter(tmy_data.DryBulb, label='Temperature (deg C)')


Out[32]:
(<matplotlib.figure.Figure at 0x10b35a20>,
 <matplotlib.axes._subplots.AxesSubplot at 0x10b37908>)

In [33]:
sapm_sd_scatter(tmy_data.DNI, label='DNI (W/m**2)')


Out[33]:
(<matplotlib.figure.Figure at 0x10e63080>,
 <matplotlib.axes._subplots.AxesSubplot at 0x10bf0a90>)

In [34]:
sapm_sd_scatter(tmy_data.AOD, label='AOD')


Out[34]:
(<matplotlib.figure.Figure at 0x10e63860>,
 <matplotlib.axes._subplots.AxesSubplot at 0x10b89a90>)

In [51]:
sapm_sd_scatter(tmy_data.Wspd, label='Wind speed', vmax=10)


Out[51]:
(<matplotlib.figure.Figure at 0x11c157390>,
 <matplotlib.axes._subplots.AxesSubplot at 0x11e88a490>)

Notice the use of the vmax keyword argument in the above example. The **kwargs pattern allows us to easily pass non-specified arguments to nested functions.


In [52]:
def sapm_other_scatter(c_data, x_data, clabel=None, xlabel=None, aspect_equal=False, **kwargs):
    """Display a scatter plot of SAPM p_ac vs. something else.
    
    You need to re-execute this cell if you re-run the p_ac calculation.
    
    Parameters
    ----------
    c_data : array-like
        Determines the color of each point on the scatter plot.
        Must be same length as p_acs.
    x_data : array-like
        
    kwargs passed to ``scatter``.
    
    Returns
    -------
    tuple of fig, ax objects
    """
    
    fig = plt.figure(figsize=(12,12))
    
    if aspect_equal:
        ax = fig.add_subplot(111, aspect='equal')
    else:
        ax = fig.add_subplot(111)
    sc = ax.scatter(x_data, p_acs['sapm'], c=c_data, alpha=1, cmap=mpl.cm.YlGnBu_r, **kwargs)  
    ax.set_xlim(0, None)
    ax.set_ylim(0, None)
    ax.set_xlabel('{}'.format(xlabel))
    ax.set_ylabel('Sandia model power (W)')
    fig.colorbar(sc, label='{}'.format(clabel), shrink=0.75)
    
    return fig, ax

In [53]:
sapm_other_scatter(tmy_data.DryBulb, tmy_data.GHI, clabel='Temperature (deg C)', xlabel='GHI (W/m**2)')


Out[53]:
(<matplotlib.figure.Figure at 0x11e891110>,
 <matplotlib.axes._subplots.AxesSubplot at 0x11e8bb0d0>)

Next, we will assume that the SAPM model is representative of the real world performance so that we can use scipy's optimization routine to derive simulated PVUSA coefficients. You will need to install scipy to run these functions.

Here's one PVUSA reference:

http://www.nrel.gov/docs/fy09osti/45376.pdf


In [54]:
def pvusa(pvusa_data, a, b, c, d):
    """
    Calculates system power according to the PVUSA equation
    
    P = I * (a + b*I + c*W + d*T)
    
    where
    P is the output power,
    I is the plane of array irradiance,
    W is the wind speed, and
    T is the temperature
    
    Parameters
    ----------
    pvusa_data : pd.DataFrame
        Must contain the columns 'I', 'W', and 'T'
    a : float
        I coefficient
    b : float
        I*I coefficient
    c : float
        I*W coefficient
    d : float
        I*T coefficient
    
    Returns
    -------
    power : pd.Series
        Power calculated using the PVUSA model.
    """
    return pvusa_data['I'] * (a + b*pvusa_data['I'] + c*pvusa_data['W'] + d*pvusa_data['T'])

In [55]:
from scipy import optimize

In [56]:
pvusa_data = pd.DataFrame()
pvusa_data['I'] = poa_irrad.poa_global
pvusa_data['W'] = tmy_data.Wspd
pvusa_data['T'] = tmy_data.DryBulb

In [57]:
popt, pcov = optimize.curve_fit(pvusa, pvusa_data.dropna(), p_acs.sapm.values, p0=(.0001,0.0001,.001,.001))
print('optimized coefs:\n{}'.format(popt))
print('covariances:\n{}'.format(pcov))


optimized coefs:
[  2.09706292e-01   5.45378922e-06   9.35100724e-04  -1.37507303e-03]
covariances:
[[  1.34031678e-07  -9.74577537e-11  -7.82618534e-09  -2.31707016e-09]
 [ -9.74577537e-11   1.87691716e-13  -7.69298228e-13  -2.50824263e-12]
 [ -7.82618534e-09  -7.69298228e-13   1.25182333e-09   1.91488592e-10]
 [ -2.31707016e-09  -2.50824263e-12   1.91488592e-10   3.76869810e-10]]

In [58]:
power_pvusa = pvusa(pvusa_data, *popt)

fig, ax = sapm_other_scatter(tmy_data.DryBulb, power_pvusa, clabel='Temperature (deg C)',
                             aspect_equal=True, xlabel='PVUSA (W)')

maxmax = max(ax.get_xlim()[1], ax.get_ylim()[1])
ax.set_ylim(None, maxmax)
ax.set_xlim(None, maxmax)
ax.plot(np.arange(maxmax), np.arange(maxmax), 'r')


Out[58]:
[<matplotlib.lines.Line2D at 0x11eb64590>]

In [ ]: