Update intro tutorial (#1116)

Make the procedural and OO results match exactly and highlight the use
of `iotools` to download meteorological data.

Author: Peque

docs/sphinx/source/introtutorial.rst

@@ -27,19 +27,18 @@ configuration at a handful of sites listed below.
 
 .. ipython:: python
 
+    import pvlib
     import pandas as pd
     import matplotlib.pyplot as plt
 
-    naive_times = pd.date_range(start='2015', end='2016', freq='1h')
-
     # very approximate
     # latitude, longitude, name, altitude, timezone
-    coordinates = [(30, -110, 'Tucson', 700, 'Etc/GMT+7'),
-                   (35, -105, 'Albuquerque', 1500, 'Etc/GMT+7'),
-                   (40, -120, 'San Francisco', 10, 'Etc/GMT+8'),
-                   (50, 10, 'Berlin', 34, 'Etc/GMT-1')]
-
-    import pvlib
+    coordinates = [
+        (30, -110, 'Tucson', 700, 'Etc/GMT+7'),
+        (35, -105, 'Albuquerque', 1500, 'Etc/GMT+7'),
+        (40, -120, 'San Francisco', 10, 'Etc/GMT+8'),
+        (50, 10, 'Berlin', 34, 'Etc/GMT-1'),
+    ]
 
     # get the module and inverter specifications from SAM
     sandia_modules = pvlib.pvsystem.retrieve_sam('SandiaMod')
@@ -48,9 +47,29 @@ configuration at a handful of sites listed below.
     inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208__208V_']
     temperature_model_parameters = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
 
-    # specify constant ambient air temp and wind for simplicity
-    temp_air = 20
-    wind_speed = 0
+
+Let's download the typical meteorological year weather data from PVGIS, which
+includes irradiation, temperature and wind speed.  Note that PVGIS uses
+different naming conventions, so it is required to rename the weather data
+variables before using them. PVGIS weather data is already UTC-localized, so
+conversion to local timezone is optional.
+
+.. ipython:: python
+
+    variables_translation = {
+        "Gb(n)": "dni",
+        "G(h)": "ghi",
+        "Gd(h)": "dhi",
+        "T2m": "temp_air",
+        "WS10m": "wind_speed",
+    }
+    tmys = []
+    for location in coordinates:
+        latitude, longitude, name, altitude, timezone = location
+        weather = pvlib.iotools.get_pvgis_tmy(latitude, longitude)[0]
+        weather = weather.rename(columns=variables_translation)
+        weather.index.name = "utc_time"
+        tmys.append(weather)
 
 
 Procedural
@@ -69,41 +88,60 @@ to accomplish our system modeling goal:
 
     energies = {}
 
-    for latitude, longitude, name, altitude, timezone in coordinates:
-        times = naive_times.tz_localize(timezone)
+    for location, weather in zip(coordinates, tmys):
+        latitude, longitude, name, altitude, timezone = location
         system['surface_tilt'] = latitude
-        solpos = pvlib.solarposition.get_solarposition(times, latitude, longitude)
-        dni_extra = pvlib.irradiance.get_extra_radiation(times)
+        solpos = pvlib.solarposition.get_solarposition(
+            time=weather.index,
+            latitude=latitude,
+            longitude=longitude,
+            altitude=altitude,
+            temperature=weather["temp_air"],
+            pressure=pvlib.atmosphere.alt2pres(altitude),
+        )
+        dni_extra = pvlib.irradiance.get_extra_radiation(weather.index)
         airmass = pvlib.atmosphere.get_relative_airmass(solpos['apparent_zenith'])
         pressure = pvlib.atmosphere.alt2pres(altitude)
         am_abs = pvlib.atmosphere.get_absolute_airmass(airmass, pressure)
-        tl = pvlib.clearsky.lookup_linke_turbidity(times, latitude, longitude)
-        cs = pvlib.clearsky.ineichen(solpos['apparent_zenith'], am_abs, tl,
-                                     dni_extra=dni_extra, altitude=altitude)
-        aoi = pvlib.irradiance.aoi(system['surface_tilt'], system['surface_azimuth'],
-                                   solpos['apparent_zenith'], solpos['azimuth'])
-        total_irrad = pvlib.irradiance.get_total_irradiance(system['surface_tilt'],
-                                                            system['surface_azimuth'],
-                                                            solpos['apparent_zenith'],
-                                                            solpos['azimuth'],
-                                                            cs['dni'], cs['ghi'], cs['dhi'],
-                                                            dni_extra=dni_extra,
-                                                            model='haydavies')
-        tcell = pvlib.temperature.sapm_cell(total_irrad['poa_global'],
-                                            temp_air, wind_speed,
-                                            **temperature_model_parameters)
+        aoi = pvlib.irradiance.aoi(
+            system['surface_tilt'],
+            system['surface_azimuth'],
+            solpos["apparent_zenith"],
+            solpos["azimuth"],
+        )
+        total_irradiance = pvlib.irradiance.get_total_irradiance(
+            system['surface_tilt'],
+            system['surface_azimuth'],
+            solpos['apparent_zenith'],
+            solpos['azimuth'],
+            weather['dni'],
+            weather['ghi'],
+            weather['dhi'],
+            dni_extra=dni_extra,
+            model='haydavies',
+        )
+        cell_temperature = pvlib.temperature.sapm_cell(
+            total_irradiance['poa_global'],
+            weather["temp_air"],
+            weather["wind_speed"],
+            **temperature_model_parameters,
+        )
         effective_irradiance = pvlib.pvsystem.sapm_effective_irradiance(
-            total_irrad['poa_direct'], total_irrad['poa_diffuse'],
-            am_abs, aoi, module)
-        dc = pvlib.pvsystem.sapm(effective_irradiance, tcell, module)
+            total_irradiance['poa_direct'],
+            total_irradiance['poa_diffuse'],
+            am_abs,
+            aoi,
+            module,
+        )
+        dc = pvlib.pvsystem.sapm(effective_irradiance, cell_temperature, module)
         ac = pvlib.inverter.sandia(dc['v_mp'], dc['p_mp'], inverter)
         annual_energy = ac.sum()
         energies[name] = annual_energy
 
     energies = pd.Series(energies)
 
     # based on the parameters specified above, these are in W*hrs
-    print(energies.round(0))
+    print(energies)
 
     energies.plot(kind='bar', rot=0)
     @savefig proc-energies.png width=6in
@@ -150,28 +188,35 @@ by examining the parameters defined for the module.
     from pvlib.location import Location
     from pvlib.modelchain import ModelChain
 
-    system = PVSystem(module_parameters=module,
-                      inverter_parameters=inverter,
-                      temperature_model_parameters=temperature_model_parameters)
+    system = PVSystem(
+        module_parameters=module,
+        inverter_parameters=inverter,
+        temperature_model_parameters=temperature_model_parameters,
+    )
 
     energies = {}
-    for latitude, longitude, name, altitude, timezone in coordinates:
-        times = naive_times.tz_localize(timezone)
-        location = Location(latitude, longitude, name=name, altitude=altitude,
-                            tz=timezone)
-        weather = location.get_clearsky(times)
-        mc = ModelChain(system, location,
-                        orientation_strategy='south_at_latitude_tilt')
-        # model results (ac, dc) and intermediates (aoi, temps, etc.)
-        # assigned as mc object attributes
-        mc.run_model(weather)
-        annual_energy = mc.results.ac.sum()
+    for location, weather in zip(coordinates, tmys):
+        latitude, longitude, name, altitude, timezone = location
+        location = Location(
+            latitude,
+            longitude,
+            name=name,
+            altitude=altitude,
+            tz=timezone,
+        )
+        mc = ModelChain(
+            system,
+            location,
+            orientation_strategy='south_at_latitude_tilt',
+        )
+        results = mc.run_model(weather)
+        annual_energy = results.ac.sum()
         energies[name] = annual_energy
 
     energies = pd.Series(energies)
 
     # based on the parameters specified above, these are in W*hrs
-    print(energies.round(0))
+    print(energies)
 
     energies.plot(kind='bar', rot=0)
     @savefig modelchain-energies.png width=6in