Solar panel analysis pt 3: Scanning for objects

Energy loss due to spreading

Image you carry a flashlight in a dark room. When you point it directly at the wall (90° angle) you get a small circle that is intensely lit. However, if you stand close to the wall and point the flashlight somewhat sideways you light a bigger surface. The same energy emitted by the flashlight is now spread out and each square centimeter in the lit area now receives less. We can calculate this effect with cos(pi/2 - angle). When we transform a 90° angle to radials we get pi/2 and when we enter this in the function we get the cosine of zero which is 1 or 100% spreading efficiency. In Belgium the Sun has a maximal elevation of 62° which has a spreading efficiency of 88%. During the winter solstice the Sun only reaches 15.3° at noon which has a spreading efficiency of only 26%. Fortunately, the panels themselves are somewhat inclined (20°) towards the south so that even in winter the angle of arrival is closer to 35° (57% efficiency). However, this inclination results in a more complex spreading effect because the panels now point in a particular direction. Both effects are estimated in the code below.

# angles in degrees to angles in radials
# the '_deg' part of the variable name shows it's in degrees
deg_to_rad <- function(angle_deg){pi * angle_deg / 180}

# vertical spreading effect
get_vertical_spreading <- function(angle_deg){
  cos(pi / 2 - deg_to_rad(angle_deg))
  }

# When the Sun points to the south we should add 20° to the
# solar elevation to get the correct angle since the panels
# are inclined by 20° towards the south.
# However if the Sun were to shine from the north we should
# detract 20° from the elevation. The negative cosine 
# should do just this for us.
adjust_horizontal_angle <- function(angle_deg, azimuth_deg){
  angle_deg + 20 * -cos(deg_to_rad(azimuth_deg))
}

Energy loss due to traveling through the atmosphere

A second, less straightforward, effect that decreases the energy transfer from the Sun to the panels at lower solar elevations is the fact that the distance traveled through the atmosphere becomes larger. As the light rays make their way through this mass of air they often collide with the molecules and become scattered. It is this effect that makes sunsets so beautiful and makes the sky above above us appear blue instead of black. Some clever people have studied this effect in great detail and we can take their equations from Wikipedia to learn that at a 90° angle the energy per square meter is about 1051 Watts. At 62° (Belgium summer solstice) we get 1019 Watts/m² which is still quite good. However at 15.3° in winter we fall back to 628 Watts/m². While we lose a lot of energy due to this scattering effect a little bit may be recuperated. Some of the scattered light from rays passing overhead now reaches our panels.

# calulate Air Mass (AM) according to Karsten & Young 
# https://en.wikipedia.org/wiki/Air_mass_(solar_energy)#math_I.1
# https://en.wikipedia.org/wiki/Solar_irradiance
get_air_mass <- function(angle_deg){
  1 / (cos(deg_to_rad(angle_deg)) + 0.506 * (96.08 - angle_deg) ^ -1.6364)
}
 
# this function calculates how much power (W/m²) 
# is left in Sun rays reaching the surface
# the only effect weakening solar power is air mass at this point
# incoming solar power on top of atmosphere = 1366 W/m²
get_solar_intensity <- function(angle_deg){
  # extract from 90 to get correct z angle (see wikipedia)
  1.1 * 1366 * 0.7 ^ get_air_mass(90 - angle_deg) ^ 0.678
}

# Raw estimation of the scattered light reaching the panels
get_scatterlight_intensity <- function(angle_deg){
  
  # The solar intensity of rays passing overhead 
  solar_intensity <- get_solar_intensity(angle_deg + 3)
  # The percentage of max energy (at 90°) that remains
  pct_intensity <- solar_intensity / get_solar_intensity(90)
  # Scattered light is a fraction of the lost energy 
  scatter_intensity <- solar_intensity * (1 - pct_intensity) * 0.15
  
  return(scatter_intensity)
}

Putting it all together

Now that we have created functions for the most important effects at play it’s time to put them together and create a master function that will estimate the maximal power yield of the solar installation when we feed it the position of the Sun.

solar_position_to_max_power <- function(angle_deg, 
                                        azimuth_angle_deg){
  
  # how much of the solar intensity is left 
  # after travelling through the air mass?
  solar_intensity <- get_solar_intensity(angle_deg)
  
  # light spreading out because of sideways spreading 
  horizontal_angle_adjusted <- adjust_horizontal_angle(angle_deg, 
                                                       azimuth_angle_deg)
  
  # panels are inclined at 20°
  vertical_spreading <- get_vertical_spreading(horizontal_angle_adjusted)
  
  surface_power <- {solar_intensity * vertical_spreading + 
                    get_scatterlight_intensity(angle_deg)}
  
  # 7 panels with 1.6 m²
  panel_surface <- 1.6 * 7 
  
  # 16% efficiency of panels
  efficieny <- 0.18
  
  # eventual power is the panel surface times 
  # the power left in the Sun's rays on the panel times
  # the efficiency of the panels
  max_plant_power <- panel_surface * surface_power * efficieny 
  
  return(max_plant_power)
}

We can now call this function on the dataset with solar positions for a full year which we created earlier on.

solar_position_max_power_df <- solar_position_df%>%
  # calculate the maximal power of the solar panels per timestamp
  # we need this to check if the maximal observed power is much lower 
  # than the maximal theoretical power
  mutate(max_power = solar_position_to_max_power(solar_elevation, 
                                                 solar_azimuth_angle))%>%
  filter(max_power >= 0)

Now let’s bring both datasets together and calculate the ratio between the power observed and the estimated maximal power.

combined_df <- solarplant_df%>%
  # add the solar plant data to the solar position data
  # standardise the time to 1970 to allow joing the solar position dataset
  mutate(timestamp_basic = lubridate::origin + (yday(timestamp) - 1) * 3600 * 24 + 
                           hour(timestamp) * 3600 + minute(timestamp) * 60)%>%
  right_join(solar_position_max_power_df, by = c("timestamp_basic" = "timestamp"))%>%
  replace_na(list(power = 0))%>%
  # efficiency is power / max power but can't be higher than 1 (100%)
  mutate(efficiency = ifelse(power / max_power > 1, 1, power / max_power),
         # some time variables we need later on
         hour = hour(timestamp) + minute(timestamp) / 60,
         week_date = round_date(timestamp_basic, unit = "week"))

To check if our theoretical model delivers acceptable results let’s visualize both observed and estimated maximal power for some sunny days in winter and summer.

combined_df%>%
  mutate(date = as_date(timestamp))%>%
  filter(date %in% as_date(c("2017-06-01", "2016-12-04")))%>%
  ggplot(aes(x = hour(timestamp) + minute(timestamp) / 60, y = max_power))+
  geom_line()+
  geom_point(aes(y = power))+
  facet_wrap(~date)+
  labs(y = 'Power (W)',
       x = 'Time of day (h)')+
  theme_minimal()+
  theme(text = element_text(colour = "grey40"),
        strip.text.x = element_text(colour = "grey40"),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "grey40", 
                                     size = rel(0.3),
                                     fill = NA))

As you can see our theoretical estimation fits the observed values quite well both in winter and summer. While not perfect, this estimation should now allow us to detect solar positions that never result in maximal efficiency values due to an obstacle blocking the light.

Visualizing maximal efficieny per solar position

Most of the time when the panels yield less than what is theoretically possible, clouds are spoiling the party. To filter out this noise and let the effect of stationary objects emerge we will look at the maximal efficiency for each time of the day, per week of the year. Since we have 2 years of solar data in the dataset we now have registered 14 days for each week of the year. For clouds to appear as a stationary object in our analysis it would have to be cloudy on a particular time of the day for all 14 days measured. Even for Belgium this is rarely the case.

max_efficiency_df <- combined_df%>%
  # for each week of the year we hope the Sun shined at 
  # least once per 10 minute interval this way we can get 
  # an accurate maximal efficiency estimate by grouping 
  # on the week and hour variables and calculating
  # the maximal values for all relevant variables
  group_by(week_date, hour)%>%
  summarise(solar_azimuth_angle = max(solar_azimuth_angle),
            solar_elevation = max(solar_elevation),
            efficiency = max(efficiency))%>%
  ungroup()

max_efficiency_df%>%
  ggplot(aes(x = solar_azimuth_angle, 
             y = solar_elevation, 
             colour = efficiency, 
             alpha = efficiency,
             size = efficiency))+
  geom_point()+
  scale_color_gradient2(low = "darkgreen", 
                        mid = "lightblue", 
                        high = "lightblue", 
                        midpoint = 0.7, 
                        labels = scales::percent)+
  scale_alpha_continuous(range = c(1, 0.2), 
                        labels = scales::percent)+
  scale_size_continuous(range = c(8, 1), 
                         labels = scales::percent)+
  scale_x_continuous(breaks = seq(45, 315, 45), 
                     expand = c(0, 0),
                     labels = c("Northeast", "East", 
                                "Southeast", "South", 
                                "Southwest", "West",
                                "Northwest"))+
  scale_y_continuous(expand = c(0, 0), limits = c(0, 65))+
  labs(y = "Solar elevation (°)")+
  theme(text = element_text(colour = "grey40"),
        panel.background = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        legend.key = element_blank(),
        legend.justification = 'top')+
  guides(colour = guide_legend(title.position = "top", 
                               reverse = T,
                               ncol = 1),
         alpha = guide_legend(title.position = "top", 
                              reverse = T,
                               ncol = 1),
         size = guide_legend(title.position = "top", 
                             reverse = T,
                               ncol = 1))

As you can see I’ve taken some artistic liberty while creating this plot. Low maximal efficiencies are colored in green and high values maximal values get a blue color + higher transparency and smaller circles.

To verify this result I got up on my roof and took a panoramic picture. Let’s first look at this picture without the data.

And now let’s use it as the background for our plot.

library(jpeg)
library(grid)
img <- jpeg::readJPEG(rt$find_file("static/figs/roof_panorama.jpg"))

g <- rasterGrob(img, interpolate = TRUE) 

max_efficiency_df%>%
  ggplot(aes(x = solar_azimuth_angle, 
             y = solar_elevation, 
             colour = efficiency, 
             alpha = efficiency, 
             size = efficiency))+
  annotation_custom(g, xmin = -50, xmax = 320, 
                    ymin = -15, ymax = 65)+
  geom_point()+
  scale_color_gradient2(low = "darkgreen", 
                        mid = "lightblue", 
                        high = "lightblue", 
                        midpoint = 0.7, 
                        labels = scales::percent)+
  scale_alpha_continuous(range = c(1, 0.2), 
                        labels = scales::percent)+
  scale_size_continuous(range = c(8, 1), 
                         labels = scales::percent)+
  scale_x_continuous(breaks = seq(45, 315, 45), 
                     labels = c("Northeast", "East", 
                                "Southeast", "South", 
                                "Southwest", "West",
                                "Northwest"))+
  scale_y_continuous(expand = c(0, 0), limits = c(0, 65))+
  labs(y = "Solar elevation (°)")+
  theme(text = element_text(colour = "grey40"),
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = 'none',
        axis.title.x = element_blank())

Et voilà! The objects we detected do indeed match the trees we suspected all along. We’ve successfully used the Sun to scan the surroundings for objects using solar panel data! The panels estimate the tree to the east to be taller but this is because the picture was taken on the side of the roof furthest from the tree. To the panels on the other side of the roof the tree seems taller.

I expect this to be the final post on this dataset and hope you enjoyed the analysis as much as I did. Time for a new pet project. Cheers!