Using shading function on 15-minute data

[matsuobasho]

It is rare for real-world commercial, non-test pv plant data to be at a 1-minute resolution. Is there a way to adopt the shading function to work with 15-minute time intervals?

[cwhanse]

@matsuobasho for this particular shading function? No. We're actually considering withdrawing the algorithm from the library. After using it for many locations, it appears this algorithm requires extended days of nearly-clear conditions to work well.

We're very interested to hear about other algorithms that might help identify structural shading.

[williamhobbs]

@cwhanse If it's of interest, I can share details on what I did a few years ago on this topic. I think you and I discussed this at the time, so there may not be anything new here. And additional processing would be needed to get a binary shade mask. Maybe conventional image processing techniques would work on my final images, especially if they had been normalized to clear sky to remove the gradient.

(from PVPMC 2017 in ABQ, https://www.slideshare.net/sandiaecis/04-final-hobbs-lave-wvm-solar-portfolios-pvpmc)

Here's a comparison of Solmetric SunEye images/charts and results using a pyranometer and power from a single module for this rooftop system, https://goo.gl/maps/MHmy6ci1cAMqZXt47:

[AdamRJensen]

@williamhobbs The irradiance you plot, is that the average for each bin? I have often found it useful to plot the maximum for each bin or perhaps 99th percentile (this basically filters out cloudy periods). When doing that the differences are much clearer/starker contrast.

Nice work by the way ☺️

[cwhanse]

@AdamRJensen what is meant by "bin"? Are you thinking to bin the irradiance (e.g., by interval of 50 W/m2) before making the heat plots? If so, what is the extent of a bin: a time interval x a few days perhaps?

[williamhobbs]

@cwhanse I put the irradiance values in m*m degree bins (azimuth and elevation) and then replaced values with the nth percentile from each bin. Sort of.

@AdamRJensen It's a percentile.

My code was poorly commented and I haven't spent the time to go back through it all the way and figure out what I did, but it looks like I was using the 80th percentile in one step. That seems low - if I did this from scratch now I would probably start with 95th-99th, so I'm not sure what was going on...

A copy of my MATLAB code is at the bottom of this comment, in case it gives anyone any ideas.

I haven't fully thought this through, but it seems like there could be advantages to processing in both time-of-day/day-of-year coordinates and and azimuth/elevation coordinates. For example:

• In az/el bins, replace each value with the nh percentile
• But the results back into time/day coordinates and apply a 1D filter to remove low values that go across time of day (horzintally) but not day of year (vertically)
• put it back in az/el format and apply morphological steps

printFigs = 0;

if exist('PV_Data_Shade_Sensor_Size_Processed.mat','file')
    load('PV_Data_Shade_Sensor_Size_Processed.mat')
else
    
    load 'PV_Data_Shade_Sensor_Size.mat'
    
    prctileThresh = .8;
    cloudThresh = .6;
    
    for i=1:length(siteNames)
        
        PV.(siteNames{i}).indDates = true(length(PV.(siteNames{i}).dt),1);%PV.(siteNames{i}).dt > startDate & PV.(siteNames{i}).dt < stopDate;
        
        PV.(siteNames{i}).irrFiltered=PV.(siteNames{i}).POA;
        PV.(siteNames{i}).isCloud = false(length(PV.(siteNames{i}).irrFiltered),1);
        avgIrr = zeros(length(PV.(siteNames{i}).irrFiltered),1);
        PV.(siteNames{i}).prcIrr = zeros(length(PV.(siteNames{i}).irrFiltered),1);
        
        PV.(siteNames{i}).minOfDay = PV.(siteNames{i}).Time.hour*60+PV.(siteNames{i}).Time.minute+1;
        PV.(siteNames{i}).dayOfPeriod = floor(PV.(siteNames{i}).dt-PV.(siteNames{i}).dt(1))+1; %day of year
        
        for j=min(PV.(siteNames{i}).dayOfPeriod(PV.(siteNames{i}).indDates)):max(PV.(siteNames{i}).dayOfPeriod(PV.(siteNames{i}).indDates))
            for k=200:1200% 1440
                               
                % For n = prctileThresh*100, calculate nth percentile of
                % irradiance for all measurements where SunAz and SunEl are +/-
                % m degrees from current day of year and minute of day. This is
                % a proxy for clear sky irradiance, leaving in fixed shading.
                m=0.5;
                currSunEl = PV.(siteNames{i}).SunEl(PV.(siteNames{i}).dayOfPeriod == j & PV.(siteNames{i}).minOfDay == k);
                
                if currSunEl >= 0 % if it is daytime
                    currSunAz = PV.(siteNames{i}).SunAz(PV.(siteNames{i}).dayOfPeriod == j & PV.(siteNames{i}).minOfDay == k);
                    PV.(siteNames{i}).indAngle = PV.(siteNames{i}).SunAz>=currSunAz-m & PV.(siteNames{i}).SunAz<=currSunAz+m & PV.(siteNames{i}).SunEl>=currSunEl-m & PV.(siteNames{i}).SunEl<=currSunEl+m;
                    PV.(siteNames{i}).indAngle = PV.(siteNames{i}).indAngle & PV.(siteNames{i}).indDates;
                    PV.(siteNames{i}).prcIrr(PV.(siteNames{i}).dayOfPeriod==j&PV.(siteNames{i}).minOfDay==k)=percentile(PV.(siteNames{i}).POA(PV.(siteNames{i}).indAngle),prctileThresh);
                end
                
            end
            
            if rem(j,5)==0
                disp({'site = ', i,', day = ',j})
            end
        end
        
        PV.(siteNames{i}).irrFiltered(PV.(siteNames{i}).POA<PV.(siteNames{i}).prcIrr*cloudThresh) = PV.(siteNames{i}).prcIrr(PV.(siteNames{i}).POA<PV.(siteNames{i}).prcIrr*cloudThresh);
        PV.(siteNames{i}).isCloud(PV.(siteNames{i}).POA<PV.(siteNames{i}).prcIrr*cloudThresh) = 1;
        
    end
    
    save PV_Data_Shade_Sensor_Size_Processed.mat -v7.3 % tables can only be saved with the -v7.3 switch
    
end

for i=1:length(siteNames)
    %% Day by minute plot
    figure;
    imagesc([0:1/1440:1],[180,364],reshape(PV.(siteNames{i}).POA,1440,[])')
    datetick('x','hh:MM')
    datetick('y')
    xlim([0 .99])
    ylim([180,364])
    h=colorbar;
    ylabel(h, 'Irradiance W/m^2)')
    ylabel('Month')
    xlabel('Time of Day (CST)')
    title(siteNames{i})
    
    %http://www.mathworks.com/matlabcentral/newsreader/view_thread/169200
    set(gcf,'PaperPositionMode','auto','units','inches','pos', [3 3 8 4])
    
    %set background to white
    set(gcf,'color','w');
    
    if printFigs
        print(['Unfiltered_Day_by_Minute_Shade_Plot_',siteNames{i},'.png'],'-dpng','-r300');
    end
    
    %% unfiltered shade plot
    figure
    scatter(PV.(siteNames{i}).SunAz(PV.(siteNames{i}).indDates),PV.(siteNames{i}).SunEl(PV.(siteNames{i}).indDates),5*ones(length(PV.(siteNames{i}).POA(PV.(siteNames{i}).indDates)),1),PV.(siteNames{i}).POA(PV.(siteNames{i}).indDates),'filled')
    ylim([0 90])
    xlabel('Azimuth (Deg)')
    ylabel('Elevation (Deg)')
    h=colorbar;
    ylabel(h,'irradiance (W/m^2)')
    title(siteNames{i})
    
    %http://www.mathworks.com/matlabcentral/newsreader/view_thread/169200
    set(gcf,'PaperPositionMode','auto','units','inches','pos', [3 3 8 4])
    
    %set background to white
    set(gcf,'color','w');
    
    if printFigs
        print(['Unfiltered_Shade_Plot_',siteNames{i},'.png'],'-dpng','-r300');
    end
    
    %% filtered shade plot
    figure
    scatter(PV.(siteNames{i}).SunAz(PV.(siteNames{i}).indDates),PV.(siteNames{i}).SunEl(PV.(siteNames{i}).indDates),5*ones(length(PV.(siteNames{i}).POA(PV.(siteNames{i}).indDates)),1),PV.(siteNames{i}).irrFiltered(PV.(siteNames{i}).indDates),'filled')
    ylim([0 90])
    xlabel('Azimuth (Deg)')
    ylabel('Elevation (Deg)')
    h=colorbar;
    ylabel(h,'irradiance (W/m^2)')
    title(siteNames{i})
    
    %http://www.mathworks.com/matlabcentral/newsreader/view_thread/169200
    set(gcf,'PaperPositionMode','auto','units','inches','pos', [3 3 8 4])
    
    %set background to white
    set(gcf,'color','w');
    
    if printFigs
        print(['Shade_Plot_',siteNames{i},'.png'],'-dpng','-r300');
    end
end

[cwhanse]

put it back in az/el format and apply morphological steps

I don't know how to form the morphological elements in a radial coordinate system, I suppose it's possible. In the Cartesian (time-of-day, day-of-year) image it's easy to define the structuring elements. The structuring elements are akin to sliding windows, and the algorithm looks for adjacent elements that meet some threshold, to assemble connected structure in the image. Maybe I'm overthinking this and the same rectangular structuring elements would work.

[williamhobbs]

Good point. I never ended up with a true "image" in radial (az, el) coordinates. I just did a scatter plot, where markers got stacked on top of each other, so plotting order could change the final appearance. A proper image transformation would be needed, I think.

For the data from the distribution pole, I made a few changes from the code I posted previously. This code might be easier to follow.

"DPV" is a MATLAB struct with data and metadata for a bunch of sites, "Hoov6" is the one used here.

load 'DPV_Data_Shade.mat'

startDate=datenum([2012 1 1]);
stopDate=datenum([2012 12 21]);
indDates = DPV.Hoov6.dt > startDate & DPV.Hoov6.dt < stopDate;
   
irr=DPV.Hoov6.POA;
irrFiltered=irr;
prctileThresh = .8;
cloudThresh = .6;
isCloud = false(length(irrFiltered),1);
avgIrr = zeros(length(irrFiltered),1);
prcIrr = zeros(length(irrFiltered),1);

minOfDay = DPV.Hoov6.Time.hour*60+DPV.Hoov6.Time.minute+1;
dayOfPeriod = floor(DPV.Hoov6.dt-DPV.Hoov6.dt(1))+1; %day of year

for j=min(dayOfPeriod(indDates)):max(dayOfPeriod(indDates))
    for k=200:1200% 1440
        
        % For n = prctileThresh*100, calculate nth percentile of
        % irradiance for all measurements where SunAz and SunEl are +/-
        % m degrees from current day of year and minute of day. This is
        % a proxy for clear sky irradiance, leaving in fixed shading.
        m=0.3;
        currSunEl = DPV.Hoov6.SunEl(dayOfPeriod == j & minOfDay == k);
        
        if currSunEl >= 0 % if it is daytime
            currSunAz = DPV.Hoov6.SunAz(dayOfPeriod == j & minOfDay == k);
            indAngle = DPV.Hoov6.SunAz>=currSunAz-m & DPV.Hoov6.SunAz<=currSunAz+m & DPV.Hoov6.SunEl>=currSunEl-m & DPV.Hoov6.SunEl<=currSunEl+m;
            indAngle = indAngle & indDates;
            prcIrr(dayOfPeriod==j&minOfDay==k)=percentile(irr(indAngle),prctileThresh);
        end
        
    end
    
    if rem(j,5)==0
        disp({'site = ', i,', day = ',j})
    end
end

irrFiltered(irr<prcIrr*cloudThresh) = prcIrr(irr<prcIrr*cloudThresh);
isCloud(irr<prcIrr*cloudThresh) = 1;

%% filtered shade plot
figure
scatter(DPV.Hoov6.SunAz(indDates),DPV.Hoov6.SunEl(indDates),5*ones(length(DPV.Hoov6.POA(indDates)),1),irrFiltered(indDates),'filled')
ylim([0 90])

ylabel('Azimuth (deg)')
xlabel('Elevation (deg)')

h=colorbar; 
ylabel(h,'irradiance (W/m^2)')

[williamhobbs]

I just saw this paper that seems relevant: https://doi.org/10.1002/pip.3654

[AdamRJensen]

The methodology for detecting shading proposed in Lorenz, E. (2022) is also rather interesting.

[cwhanse]

Two different shading sources: near-field from objects and horizon. The horizon shading in Lorenz et al. requires identifying clear-sky conditions (they use satellite-derived irradiance); is it reasonable to expect users to obtain that? It's not a blocker for implementing that horizon shade detector (it's a nice algorithm with relatively simple steps, the machine learning is just clustering), but would need to be explicit about the required inputs.

[williamhobbs]

@cwhanse, if other parts of pvanalytics start making use of "external" data (e.g., #143), then I could see a workflow that uses satellite data (or ERA5, MERRA2) to detect clear-sky being reasonable.

[cwhanse]

We had imagined that pvanalytics would also contain workflows that assemble base-layer functions to accomplish some purpose. Workflows would be the analogue of the ModelChain object in pvlib-python.