Renewable Energy Analysis in R

This tutorial explores analysis and integration of renewable energy data from a variety of sources in R.

Energy Demand

The US Energy Information Agency (EIA) was created by the 1974 Federal Energy Administration (FEA) Act (P.L. 93-275, 15 USC 761) as a response to the early 1970s "energy crisis" and the need for additional Federal initiatives to collect, evaluate, analyze, and disseminate energy-related information (EIA 2021).

State Energy Profiles

EIA's State Energy Data System (SEDS) is a compilation of time series data on state energy production, consumption, prices, and expenditures. The objective is to provide data that is defined as consistently as possible over time and across sectors for analysis and forecasting purposes (EiA 2021).

A GeoJSON file of selected SEDS variables is available here.


states = st_read("2019-state-energy.geojson")


	 [1] "ST"                              "Name"                           
	 [3] "GEOID"                           "AFFGEOID"                       
	 [5] "Square.Miles.Land"               "Square.Miles.Water"             
	 [7] "State.Name"                      "Population.MM"                  
	 [9] "Civilian.Labor.Force.MM"         "GDP.B"                          
	[11] "GDP.Per.Capita"                  "GDP.Manufacturing.MM"           
	[13] "Personal.Income.Per.Capita"      "VMT.MM"                         
	[15] "Farm.Land.MM.Acres"              "Average.Temperature"            
	[17] "Precipitation.Annual"            "Consumption.Coal.B.BTU"         
	[19] "Consumption.Jet.Fuel.B.BTU"      "Consumption.Gasoline.B.BTU"     
	[21] "Consumption.Total.B.BTU"         "Consumption.Per.Capita.MM.BTU"  
	[23] "Production.Coal.B.BTU"           "Production.Gas.B.BTU"           
	[25] "Production.Oil.B.BTU"            "Production.Renewable.B.BTU"     
	[27] "Production.Total.B.BTU"          "Electricity.Interstate.Flow.GWH"
	[29] "Electricity.Average.Cost"        "Electricity.Consumption.GWH"    
	[31] "Electricity.Per.Capita.KWH"      "Electricity.Hydro.GWH"          
	[33] "Electricity.Nuclear.GWH"         "Electricity.Solar"              
	[35] "Electricity.Wind.GWH"            "CO2.Total.MM.Tonnes"            
	[37] "CO2.Per.Capita.Tonnes"           "Renewable.Standard.Type"        
	[39] "Renewable.Standard.Name"         "Renewable.Standard.Year"        
	[41] "Senators.Party"                  "geometry"                       

plot(states["Consumption.Per.Capita.MM.BTU"], breaks="quantile",
	pal=colorRampPalette(c("navy", "white", "red")))
Per capita energy consumption by state in 2019 (EIA 2021)

Energy in Illinois

In 2019, the US consumed around 100.3 quadrillion BTU (quads) of total primary energy, or around 305.4 million BTU per person. Illinois consumed around 3.96 quadrillion BTU of total primary energy, or around 312.5 million BTU per person (EIA 2021a; EIA 2021b).

illinois = states[states$ST == "IL",]

illinois = st_transform(illinois, 4326)


	bbox:           xmin: -91.51308 ymin: 36.9703 xmax: -87.49649 ymax: 42.50848
	epsg (SRID):    4269
	proj4string:    +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
	   ST     Name GEOID    AFFGEOID Square.Miles.Land Square.Miles.Water
	15 IL Illinois    17 0400000US17          55512.92           2400.229
	   State.Name Population.MM Civilian.Labor.Force.MM GDP.B GDP.Per.Capita
	15   Illinois          12.6                     6.1 897.1          71198
	   GDP.Manufacturing.MM Personal.Income.Per.Capita VMT.MM Farm.Land.MM.Acres
	15               108554                      62977 107525                 27
	   Average.Temperature Precipitation.Annual Consumption.Coal.B.BTU
	15                53.3                 41.9                 591909
	   Consumption.Jet.Fuel.B.BTU Consumption.Gasoline.B.BTU
	15                     189044                     557704
	   Consumption.Total.B.BTU Consumption.Per.Capita.MM.BTU Production.Coal.B.BTU
	15                 3958520                         312.5               1095887
	   Production.Gas.B.BTU Production.Oil.B.BTU Production.Renewable.B.BTU
	15                 2549                48045                     382793
	   Production.Total.B.BTU Electricity.Interstate.Flow.GWH
	15                2554925                         -372027
	   Electricity.Average.Cost Electricity.Consumption.GWH
	15                    28.25                      138319
	   Electricity.Per.Capita.KWH Electricity.Hydro.GWH Electricity.Nuclear.GWH
	15                      10920                   124                   98735
	   Electricity.Solar Electricity.Wind.GWH CO2.Total.MM.Tonnes
	15               251                14460               210.4
	   CO2.Per.Capita.Tonnes Renewable.Standard.Type      Renewable.Standard.Name
	15              16.69841               Mandatory Renewable Portfolio Standard
	   Renewable.Standard.Year Senators.Party                       geometry
	15                    2001     Democratic MULTIPOLYGON (((-91.51297 4...

Heat Rate

Because most current renewable energy technologies produce electricity rather than heat, in order to consider energy supply and demand in an electrified post-fossil-fuel future, we need a conversion factor to know how much electricity will be needed to satisfy future energy demand.

Assuming 100% efficiency, there are 3.412 BTU in each watt-hour of electricity. Dividing 2019 total primary energy consumption of 3,968,520 billion BTU by that conversion factor gives us 116,000 gigawatt-hours of electricity.

illinois$demand.gwh = illinois$Consumption.Total.B.BTU / 3.412

	[1] 1160176

Dividing by current electricity consumption (supplied largely by nuclear plants), we see that the renewable future will need an eight-fold increase in electrical generation.

illinois$demand.percent = 100 * illinois$demand.gwh / illinois$Electricity.Consumption.GWH

	[1] 838.7682

However, because fossil fuel plants generally only convert about one-third of the heat energy in fuel to electrical energy, EIA uses a heat rate of 9.104 BTU per watt-hour (Davis et al. 2019 C-4). This same heat rate is used with renewables (solar and wind) to keep the numbers comparable across generating sources.

Using this conversion factor gives us a somewhat less daunting need to increase electrical generation by a factor of three.

illinois$demand.gwh = illinois$Consumption.Total.B.BTU / 9.104


	[1] 434811.1

illinois$demand.percent = 100 * illinois$demand.gwh / illinois$Electricity.Consumption.GWH


	[1] 314.3538

Note that there are a number of gaps in this assessment:

However, this kind of analysis gives us some general figures that can be used to assess the magnitude of the transformation and investment needed to move to a post-fossil-fuel world.

US Census Data

This tutorial will use US Census Bureau TIGER county area polygons and demographic data from the US Census Bureau's 2015-2019 American Community Survey 5-year estimates.

A GeoJSON file of that data is available here.

counties = st_read("2015-2019-acs-counties.geojson", quiet=T)

counties = counties[counties$ST == "IL",]

plot(counties["Median.Household.Income"], breaks="quantile",
	pal=colorRampPalette(c("red", "white", "navy")))
Median household income by county in Illinois from the 2015-2019 ACS five-year estimates(USCB)

Wind Energy Potential

The Solar Energy Research, Development and Demonstration Act of 1974 established the Solar Energy Research Institute, which opened in 1977. In 1991, the institute was designated a national laboratory of the US Department of Energy and renamed to the National Renewable Energy Laboratory (NREL 2021).

NREL mission is to advance "the science and engineering of energy efficiency, sustainable transportation, and renewable power technologies and provides the knowledge to integrate and optimize energy systems" (NREL 2021). As part of their mission, NREL makes wide variety of data from their research activities available to the public.

Land-Based Wind Energy Supply Curves

NREL supply curve data is the result of assessment of renewable energy potential in terms of geographic, technical, and economic potential. Each aspect of this analysis adds a layer of complexity that is informed by related data and analysis. Factors include (NREL 2021):

Supply curve data is made available in three collections that reflect different levels of that analysis:

.zip files containing .csv files of supply curves can be downloaded from the NREL website.

NREL wind supply curves download page

Load the Supply Curve Data

The wind supply curve data is provided as a regular grid of points that can be loaded with the read.csv() function and converted to point features using the st_as_sf() function.

csv = read.csv("Reference_Access_Siting_Regime_ATB_Mid_Turbine.csv",


	[1] "X"                           "latitude"                   
	[3] "longitude"                   "area_sq_km"                 
	[5] "capacity_mw"                 "generation_mwh"             
	[7] "capacity_factor"             "wind_speed_120meters"       
	[9] "distance_to_transmission_km"

wind = st_as_sf(csv, coords=c("longitude", "latitude"), crs=4326)

plot(wind["capacity_factor"], breaks="quantile", pal=colorRampPalette(c("red", "white", "blue")))
Wind supply curve generation potential grid in mWh

Illinois Wind Supply Curve

The grid points within Illinois can be subset using the st_intersection() function.

wind = st_intersection(wind, st_geometry(illinois))

plot(wind["generation_mwh"], breaks="quantile", pal=colorRampPalette(c("red", "white", "blue")))
Wind supply curve generation potential grid for Illinois

The energy and power variables can be summed for each county by performing a spatial join using the aggregate() function.

counties$wind.capacity.mwh = aggregate(wind["generation_mwh"], counties, FUN=sum)$generation_mwh

counties$ = aggregate(wind["capacity_mw"], counties, FUN=sum)$capacity_mw

plot(counties["wind.capacity.mwh"], breaks="quantile", pal=colorRampPalette(c("red", "white", "blue")))

counties = counties[order(counties$wind.capacity.mwh, decreasing=T),]

head(st_drop_geometry(counties[,c("Name", "wind.capacity.mwh")]))

		  Name wind.capacity.mwh
	633   Iroquois          13306549
	652     McLean          12965306
	648 Livingston          12877701
	647        Lee          11467672
	605  Champaign          10743495
	645    LaSalle          10549968

Wind supply curve generation potential by Illinois county

Wind Capacity Factor

Capacity factor is the percentage of a wind turbine's nameplate capacity that it actually generates over time. Capacity factor is a function both of wind strength and consistency. Because there are few places where the wind blows strongly at all times, there will be periods where turbines are generating little or no electricity. Wind turbines are less economically feasable the lower the capacity factor.

counties$wind.capacity.factor = aggregate(wind["capacity_factor"], counties, FUN=mean)$capacity_factor

plot(counties["wind.capacity.factor"], breaks="quantile", pal=colorRampPalette(c("red", "white", "blue")))
Wind supply curve capacity factor by Illinois county

Aggregate Wind Potential

Summing the potential wind megawatt-hours, converting to BTU, and then dividing by current state consumption indicates that wind power has the potential to supply all Illinois power needs.

illinois$ = sum(counties$wind.capacity.mwh, na.rm=T)


	[1] 528758388

illinois$ = illinois$ * 9.104 / 1000


	[1] 4813816

illinois$ = 100 * illinois$ / 


	[1] 121.6065

County Wind Potential

However, the wind power resource is unevenly distributed across the state relative to population. The heavily populated counties around St. Louis and Chicago have the most notable deficits in percent of demand that can be supplied with locally-generated wind power. This would necessitate construction of long-distance power lines to bring the wind power from sparsely-populated rural areas to the population centers.

counties$wind.capacity.b.btu = counties$wind.capacity.mwh * 9.104 / 1000

counties$demand.b.btu = counties$Total.Population * 
	illinois$Consumption.Total.B.BTU / 
	sum(counties$Total.Population, na.rm=T)

counties$wind.demand.percent = 100 * counties$wind.capacity.b.btu / counties$demand.b.btu

plot(counties["wind.demand.percent"], breaks="quantile", 
	pal=colorRampPalette(c("red", "white", "blue")))
Percent of county energy demand that can be supplied by wind

Ordering by percent makes it possible to list the top and bottom counties by local wind potential.

head(st_drop_geometry(counties[,c("Name", "wind.demand.percent")]))

		 Name wind.demand.percent
	631 Henderson            2339.703
	680  Schuyler            2261.218
	671      Pope            2197.794
	600     Brown            2137.698
	683     Stark            2038.666
	628  Hamilton            1901.108
 counties = counties[order(counties$wind.demand.percent, decreasing=F),]
 head(st_drop_geometry(counties[,c("Name", "wind.demand.percent")]))

	       Name wind.demand.percent
	611    Cook          0.07601945
	617  DuPage          0.20187488
	644    Lake          1.22054894
	694    Will         15.21790866
	640    Kane         16.77020079
	651 McHenry         28.64191445

Turbine Counts

Although offshore wind turbines can have nameplate capacity of as much as six megawatts, the largest onshore turbines are 4 mW apiece, and the 2016 average was 2 mW apiece. (EIA 2017)

Dividing the maximum installed nameplate capacity by that average size gives us the number of turbines. The average number of turbines per county is 700, although the histogram shows that some counties would need over twice that number.

counties$wind.turbines = counties$ / 2

illinois$wind.turbines = sum(counties$wind.turbines, na.rm=T)


	[1] 71581.94

print(mean(counties$wind.turbines, na.rm=T))

	[1] 701.7837

hist(counties$wind.turbines, breaks=20, 
	las=1, fg="white", col="navy", border=NA)

grid(nx=NA, ny=NULL, lty=1, col="#00000080")
Turbines per county under the maximum generation scenario

Construction Cost

The 2019 average construction cost for wind turbines was $1,391 per kWh (EIA 2021).

Multiplying that value times the number of turbines gives the estimated construction cost of $200 billion. A map of those expenditures per capita shows a significant amount of cash infusion into many rural counties. However, much of that cost would presumably be distributed to turbine manufacturers in other states or countries, and to specialized construction workers who would only stay in the county for a few years during construction. = 1391 * 1000

counties$wind.cost = counties$ *

illinois$wind.cost = sum(counties$wind.cost, na.rm=T)


	[1] 1.99141e+11

illinois$wind.cost.per.capita = illinois$wind.cost / (illinois$Population.MM * 1000000)


	[1] 15776.84

counties$wind.cost.per.capita = counties$wind.cost / counties$Total.Population

plot(counties["wind.cost.per.capita"], breaks="quantile", 
	pal=colorRampPalette(c("red", "white", "blue")))

counties = counties[order(counties$wind.cost.per.capita, decreasing=T),]

head(st_drop_geometry(counties[,c("Name", "wind.cost.per.capita", "Total.Population")]))

		 Name wind.cost.per.capita Total.Population
	671      Pope             319467.7             4203
	680  Schuyler             292870.1             6953
	631 Henderson             289115.5             6809
	600     Brown             274755.4             6628
	628  Hamilton             258636.2             8176
	683     Stark             250838.3             5447
Wind construction cost per capita

Operations and Maintenance Costs

In 2016 operations and maintenance costs for a wind farm was around $48,000 per mW (Reuters 2017).

Multiplying that times the wind capacity and dividing by population gives an estimate of per capita long-term local economic impact in each county.

counties$ = counties$ * 48000 

counties$wind.impact = counties$ / counties$Total.Population

illinois$ = sum(counties$, na.rm=T)


	[1] 6859691202

illinois$ = illinois$ / (illinois$Population.MM * 1000000)


	[1] 544.4199

illinois$median.household.income = median(counties$Median.Household.Income, na.rm=T)

plot(counties$Median.Household.Income, counties$wind.impact, fg="white", bty="n")

grid(nx=NA, ny=NULL, lty=1)

abline(v=illinois$median.household.income, col="navy")

counties = counties[order(counties$wind.impact, decreasing=T),]

head(st_drop_geometry(counties[,c("Name", "wind.impact", "Total.Population")]))

		 Name wind.impact Total.Population
	671      Pope   11024.048             4203
	680  Schuyler   10106.228             6953
	631 Henderson    9976.669             6809
	600     Brown    9481.136             6628
	628  Hamilton    8924.900             8176
	683     Stark    8655.814             5447
Long-term economic impact vs median household income by county

Operations and Maintenance Jobs

95 to 98% of jobs associated with wind farms are in the construction phase. However a significant number of technicians are needed to operate and maintain wind farms once construction is complete, although the exact number of direct and indirect jobs is highly variable across projects and countries. For the purposes of this analysis, we will use a Spanish figure of 0.2 direct operations and maintenence jobs per megawatt of installed capacity (Aldieri et al. 2019).

An estimated 28,600 maintainance and operations jobs would be created by a full deployment of wind turbines. However, plotting those jobs against long-term county unemployment rates indicates that those jobs would not necessarily be in the areas that most need them.

counties$ = round(counties$ * 0.2)

illinois$ = sum(counties$, na.rm=T)


	[1] 28633

plot(counties$Percent.Unemployment, counties$,
	las=1, fg="white", bty="n")

grid(nx=NA, ny=NULL, lty=1)

counties$Working.Age.Population = round(counties$Total.Population * 
	(100 - counties$Percent.Under.18 - counties$Percent.65.Plus) / 100)

counties = counties[order(counties$, decreasing=T),]

head(st_drop_geometry(counties[,c("Name", "", "Working.Age.Population")]))

		  Name Working.Age.Population
	633   Iroquois          715                  15853
	648 Livingston          672                  21480
	652     McLean          671                 113384
	647        Lee          593                  21080
	605  Champaign          575                 144216
	645    LaSalle          565                  65732
Wind farm operations and maintenance jobs vs. unemployment

Net Present Value

Net present value (NPV) is a method for assessing the value of an investment based on the difference between expenditures and revenues over time. However, NPV is more than simple subtraction since it takes into account the time value of money over the lifetime of a capital investment. NPV can then be used to compare a potential investment with alternative investments. If the NPV of a project is negative, then the project should be avoided because it will lose money (Investopedia 2021).

For this calculation we assume the following:

discount.rate = 0.08 = 25

wholesale.price = 50

counties$wind.revenue = (counties$wind.capacity.mwh * wholesale.price) - counties$

counties$ = sapply(1:nrow(counties), function(index)
	((-counties$wind.cost[index] / (1 + discount.rate)) +
		+ sum(sapply(, function(year)
			(counties$wind.revenue[index] / ((1 + discount.rate)^year))))))

limit = max(abs(range(counties$, na.rm=T)))

plot(counties[""], breaks = seq(-limit, limit, limit / 4),
	pal=colorRampPalette(c("red", "white", "blue")))

mean(counties$, na.rm=T)
Wind farm net present value at $50 / mWh

Installed Wind Capacity

The U.S. Wind Turbine Database

The EIA's United States Wind Turbine Database (USWTDB) provides the locations of land-based and offshore wind turbines in the United States, corresponding wind project information, and turbine technical specifications.

You can download the most recent update of the USWTDB from Download button on the EIA's US Energy Atlas page.

Downloading the USWTDB

turbines = st_read("uswtdb_v4_1_20210721.geojson")

turbines = turbines[turbines$t_state == "IL",]


	 [1] "case_id"    "t_state"    "p_name"     "p_year"     "p_tnum"    
	 [6] "p_cap"      "t_manu"     "t_model"    "t_cap"      "t_hh"      
	[11] "t_rd"       "t_rsa"      "t_ttlh"     "retrofit"   "retrofit_y"
	[16] "t_conf_atr" "t_conf_loc" "xlong"      "ylat"       "geometry"  

plot(turbines["t_cap"], pal=colorRampPalette(c("navy", "red")), reset=F)

plot(st_geometry(counties), add=T)
Wind turbines active in Illinois in 2021 by turbine size

Current Percent of Potential

Given the total installed capacity from the USWTB and the potential capacity from the supply curve, we can determine the percent of potential capacity currently being tapped.

y = aggregate(turbines, counties, FUN=length)

counties$wind.installed.turbines = ifelse(is.finite(y$t_cap), y$t_cap, 0)

y = aggregate(turbines["t_cap"], counties, FUN=sum)

counties$ = ifelse(is.finite(y$t_cap), y$t_cap / 1000, 0)

counties$wind.installed.percent = 100 * counties$ / counties$

plot(counties["wind.installed.percent"], pal=colorRampPalette(c("navy", "red")), reset=F)
Percent of potential wind capacity built by percent in Illinois counties in 2021

We can get descriptive statistics of installation from the data.

illinois$ = sum(counties$, na.rm=T)


	[1] 143163.9

illinois$wind.installed.percent = 100 * sum(counties$, na.rm=T) / 
	sum(counties$, na.rm=T)


	[1] 3.521048

counties = counties[order(counties$wind.installed.percent, decreasing=T),]

head(st_drop_geometry(counties[,c("Name", "wind.installed.percent")]))

	       Name wind.installed.percent
	652  McLean               29.90932
	616 Douglas               19.36752
	647     Lee               18.97917
	632   Henry               18.47642
	653   Macon               18.28073
	645 LaSalle               18.13521

Solar Energy Potential

Solar Energy Supply Curves

Aside from wind and hydropower, solar power is the other major renewable energy source with an extensive body of proven technology, robust ongoing research, and immense potential.

As with wind, NREL provides gridded solar supply curve data which can be downloaded from their website..

csv = read.csv("pv_reference_2020.csv",


	[1] "sc_gid"                       "capacity_factor"             
	[3] "global_horizontal_irradiance" "capacity_mw"                 
	[5] "area_sq_km"                   "latitude"                    
	[7] "longitude"                    "distance_to_transmission_km" 

solar = st_as_sf(csv, coords=c("longitude", "latitude"), crs=4326)

plot(solar["capacity_factor"], breaks="quantile", pal=colorRampPalette(c("blue", "white", "red")))
Solar supply curve capacity factor

Illinois Solar Energy Supply Curves

Subsetting to the state of Illinois, we see an expected distribtion in capacity factor, although the variance is small.

solar = st_intersection(solar, st_geometry(illinois))

plot(solar["capacity_factor"], breaks="quantile", pal=colorRampPalette(c("blue", "white", "red")))
Solar supply curve capacity factor in Illinois

The supply curve data does not include a potential annual generation figure, although one can be calculated from the capacity factor and potential capacity in mW.

solar$generation_mwh = solar$capacity_mw * solar$capacity_factor * 24 * 365

plot(solar["generation_mwh"], breaks="quantile", pal=colorRampPalette(c("blue", "white", "red")))
Solar supply generation potential in mW

County Solar Energy Supply Curves

The energy and power variables can be summed for each county by performing a spatial join using the aggregate() function.

counties$solar.capacity.mwh = aggregate(solar["generation_mwh"], counties, FUN=sum)$generation_mwh

counties$ = aggregate(solar["capacity_mw"], counties, FUN=sum)$capacity_mw

plot(counties["solar.capacity.mwh"], breaks="quantile", 
	pal=colorRampPalette(c("blue", "white", "red")))

counties = counties[order(counties$solar.capacity.mwh, decreasing=T),]

head(st_drop_geometry(counties[,c("Name", "solar.capacity.mwh")]))

	       Name solar.capacity.mwh
	652  McLean           51501691
	616 Douglas           24894366
	647     Lee           40256446
	632   Henry           45069415
	653   Macon           36563448
	645 LaSalle           44270513

Solar supply curve generation potential in annual mWh by Illinois county

Solar Capacity Factor

Capacity factor is the percentage of a solar panel's nameplate capacity that it actually generates over time. Capacity factor is a function both of solar irradience and cloud cover.

counties$solar.capacity.factor = aggregate(solar["capacity_factor"], counties, FUN=mean)$capacity_factor

plot(counties["solar.capacity.factor"], breaks="jenks", pal=colorRampPalette(c("blue", "white", "red")))
Solar capacity factor by Illinois county

Aggregate Solar Potential

Summing the potential solar megawatt-hours, converting to BTU, and then dividing by current state consumption indicates that solar power has the potential to supply all Illinois power needs.

illinois$solar.capacity.factor = mean(solar$capacity_factor, na.rm=T)


	[1] 0.1772521

illinois$solar.capacity.mwh = sum(counties$solar.capacity.mwh, na.rm=T)


	[1] 2834465381

illinois$solar.capacity.b.btu = illinois$solar.capacity.mwh * 9.104 / 1000


	[1] 25804973

illinois$solar.capacity.percent = 100 * illinois$solar.capacity.b.btu / 


	[1] 651.8844

Rooftop Installations

Unlike wind turbines that can be placed above farm fields or pasture with limited effect on use of the land for other purposes, solar panels require surface area that blocks the sun and is much more intrusive.

One common location for solar panel installation is on the roofs of buildings, especially single family homes. While residential rooftop installations are often considered unattractive and can run afoul of homeowners associations, they are nonetheless a source of surface area that does not directly interfere with the use of the land for residential purposes.

Estimation of potential surface area for solar panel installation is an area of extensive past and ongoing research. However, most methodologies require detailed building data and / or interpretation of remotely sensed data (Groppi et al. 2018).

As a rough estimate, we will use a simple model based on the following assumptions:

counties$solar.rooftops = counties$Total.Households * counties$Percent.Homeowners / 100 * 0.5

illinois$solar.rooftops = sum(counties$solar.rooftops, na.rm=T)


	[1] 1601280

counties$solar.rooftops.cost = counties$solar.rooftops * 24400

illinois$solar.rooftops.cost = sum(counties$solar.rooftops.cost, na.rm=T)


	[1] 39071236319

counties$ = counties$solar.rooftops * 8.2 / 1000

illinois$ = sum(counties$, na.rm=T)


	[1] 13130.5

counties$solar.rooftops.mwh = counties$ * counties$solar.capacity.factor * 365 * 24

counties$solar.rooftops.b.btu = counties$solar.rooftops.mwh * 9.104 / 1000

illinois$solar.rooftops.b.btu = sum(counties$solar.rooftops.b.btu, na.rm=T)


	[1] 182126

Solar Rooftop Demand Percentage

When restricting solar installation to rooftops, the percentage of demand that can be satisfied by solar drops significantly from the generous estimates in the NREL solar supply curve data.

counties$solar.rooftops.percent = 100 * counties$solar.rooftops.b.btu / counties$demand.b.btu

plot(counties["solar.rooftops.percent"], breaks="quantile", pal=colorRampPalette(c("blue", "white", "red")))

illinois$solar.rooftops.percent = 100 * sum(counties$solar.rooftops.b.btu, na.rm=T) / 
		sum(counties$demand.b.btu, na.rm=T)

	[1] 4.600692
Percent of demand that can be satisfied with rooftop solar

Net Present Value

As with wind, we can use net present value (NPV) to evaluate the cost effectiveness of solar across Illinois.

For this calculation we assume the following:

discount.rate = 0.08 = 25

wholesale.price = 190

counties$solar.rooftops.revenue = counties$solar.rooftops.mwh * wholesale.price

counties$ = sapply(1:nrow(counties), function(index)
	((-counties$solar.rooftops.cost[index] / (1 + discount.rate)) +
		+ sum(sapply(, function(year)
			(counties$solar.rooftops.revenue[index] / ((1 + discount.rate)^year))))))

limit = max(abs(range(counties$, na.rm=T)))

plot(counties[""], breaks = seq(-limit, limit, limit / 4),
	pal=colorRampPalette(c("red", "white", "blue")))

mean(counties$, na.rm=T)

	[1] 8593710
Rooftop solar net present value at $34 / mWh


Intermittency with renewable energy sources means that they are not always available. While all generating facilities need occasional downtime for maintenance, solar and wind have significant daily and seasonal periods where they generate little or no power because the sun is not shining or the wind is not blowing.

However, energy demand in an industrialized society involves a constant base load with additional cycles and surges that rise above base load in response to environmental or social conditions (such as hot summer afternoons).

There are numerous ways to address intermittency:

Energy Cycles

In order get get some sense of the amount of storage that might be needed for an Illinois powered entirely by renewables, we need fine-grained information on resource availability so we can find the peak amount of power that must be stored to provide 100% reliability to energy users. is an unfortunately named web app created by energy researchers Iain Staffell and Stefan Pfenninger that allows you to find daily and monthly capacity factors for both solar and wind at specific locations on the surface of the earth.

The source data for Staffella and Pfenninger's models is The Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2), from the Global Modeling and Assimilation Office at NASA's Goddard Space Flight Center. This dataset is a global atmospheric reanalysis that uses satellite observations since the 1980s to create a time series of grids of atmospheric and land surface conditions (NCAR 2021; Gelaro et al. 2017). This data is interpreted by the Staffell and Pfenninger's web app to estimate hourly wind and solar availability across an entire year (Pfenninger and Staffell 2016; Staffell and Pfenninger 2016).

The video below demonstrates how to download data from

Downloading data

Hourly Data

Load the wind time series into a data frame.

csv = read.csv("ninja_wind_40.1347_-88.1976_uncorrected.csv",, skip=3)


	[1] "time"        "local_time"  "electricity" "wind_speed" 

storage = csv[,c("local_time", "electricity")]

names(storage) = c("local_time", "wind")

csv = read.csv("ninja_pv_40.1347_-88.1976_uncorrected.csv",, skip=3)


	[1] "time"               "local_time"         "electricity"       
	[4] "irradiance_direct"  "irradiance_diffuse" "temperature"       

storage$solar = csv$electricity

Time Series Decomposition

Using the xts() function to create a time series and using stl() function to decompose the wind time series into seasonal, trend, and random components shows a noticeable dip in capacity factor during the summer.


storage.series = xts(storage[,c("wind", "solar")], as.POSIXlt(storage$local_time))

attr(storage.series, "frequency") = 24

decomposition = stl(as.ts(storage.series$wind), s.window=24, t.window=2400)

Time series decomposition of hourly wind capacity factor

In the opposite of wind energy, the time series decomposition of solar shows a noticeable peak during the summer.

decomposition = stl(as.ts(storage.series$solar), s.window=24, t.window=2400)

Time series decomposition of hourly solar capacity factor


One mitigation measure currently being actively developed and deployed is battery storage. Batteries are charged when there is surplus supply available, and those batteries then fill the gap when demand exceeds supply.

The primary challenge with batteries as of this writing is their expense. In 2018, utility-scale battery storage averaged around $625 per kilowatt hour (EIA 2020). Although costs can run as low as $175 per kWh, Ziegler et al. (2019) estimated costs would need to drop to $20 per kWh to be competitive with 2019 power sources.

In the simulation below, we assume a demand base load that evenly distributes total energy use across the year 24/7/365.

The bar plot shows the balance of supply and demand with blue bars on top representing surplus generation and red bars on the bottom representing unmet base load demand.

storage$ = illinois$demand.gwh * 1000 / 365 / 24

storage$ = illinois$ * storage$solar

storage$ = illinois$ * storage$wind

storage$balance = storage$ + storage$ - storage$

barplot(storage$balance, col=ifelse(storage$balance > 0, "navy", "darkred"), border=NA, las=1, space=0, tck=0)

grid(nx=NA, ny=NULL, lty=1)
Balance of renewable supply and base load demand by hour across the year

Iterative Storage Model

Adding storage to the simulation, the for() loop in this model proceeds iteratively through each hour of the year.

The overall battery capacity (in mWh) below is chosen to balance cost and unmet demand. You can modify the variable to test other possibilities.

The bar plot shows the level of battery charge on the top blue bars, and unmet demand on the bottom red bars.

storage$level = 0

storage$deficit = 0

battery.capacity = 1000000

for (x in 2:nrow(storage)) {
	if (storage$balance[x] >= 0) {
		storage$level[x] = min(storage$level[x - 1] + storage$balance[x], battery.capacity)
		storage$deficit[x] = 0
	} else if (storage$level[x - 1] >= -storage$balance[x]) {
		storage$level[x] = storage$level[x - 1] + storage$balance[x]
		storage$deficit[x] = 0 
	} else {
		storage$deficit[x] = storage$balance[x] - storage$level[x - 1]
		storage$level[x] = 0

barplot(t(storage[,c("level", "deficit")]), col=c("navy", "darkred"), border=NA, las=1, space=0, tck=0)

grid(nx=NA, ny=NULL, lty=1)
Modeled battery level vs. supply deficit

Given the wide range of costs per mWh of storage given above, we can get a range of the cost for the the battery storage modeled above.

illinois$battery.cost = list(battery.capacity * c(20, 175, 625))


[1] 2.00e+07 1.75e+08 6.25e+08

Saving Your Analysis

The three core analytical objects can be saved as .rds files so you can use them again later without having to completely re-run your analysis.

saveRDS(illinois, "2021-illinois-state-analysis.rds")
saveRDS(counties, "2021-illinois-county-analysis.rds")
saveRDS(storage, "2021-illinois-storage-analysis.rds")

They can be reloaded with the readRDS() function.

illinois = readRDS("2021-illinois-state-analysis.rds")
counties = readRDS("2021-illinois-county-analysis.rds")
storage = readRDS("2021-illinois-storage-analysis.rds")