CDS 302 - Scientific Data and Databases
Final Project (in progress): Correlation of Luminosity of Solar data with Total Solar Irradiance
Goals:
Problems with the data:
1.Download SORCE data
Below is the code I used to download the TSI data.
url = 'http://lasp.colorado.edu/data/sorce/tsi_data/daily/sorce_tsi_L3_c24h_latest.txt';
tline = 'irradiance_data.txt';
if ~exist(tline,'file') % If file is not on disk, download it.
urlwrite(url,tline);
end
fid = fopen(tline,'r');
k = 1;
while 1
tline2 = fgetl(fid);
if ~ischar(tline2),break,end;
if (length(tline2) >= 4 & tline2(1:2) == '20') %Get all the data
data(k,:) = str2num(tline2);
k = k+1;
end
end
fclose(fid);
dno = datenum(2003,2,25); %Date data began to be taken in this file
dates = [dno:dno+length(data)-1]; %Create dates vector
data(data(:,5)==0,5) = NaN; %Replace any zeroes in data with NaN.
data(:,1) = floor(data(:,1)); %Trim off decimal in date
figure(1)
plot(dates,data(:,5)) % Note that missing values (NaN) are not displayed.
hold on;
grid on;
axis tight
title('Plot of SORCE TSI data, 2012-2014','FontSize',15)
xlabel('Time (years)')
ylabel('Irradiance (W/m^2)')
datetick('x','yyyy')
xlim([734869 735964]) %Constrain x-axis to size
hold off;
- Download TSI data from SORCE
- Download years of solar videos from the Solar Dynamics Observatory directory
- Compile videos into daily average pictures
- Determine which region of the sun (say the polar regions) has luminance values that best predict the total solar irradiance (TSI)
Problems with the data:
- Variance in number of frames per video
- Black regions in some videos, eclipses
- Missing files / different datatypes
1.Download SORCE data
Below is the code I used to download the TSI data.
url = 'http://lasp.colorado.edu/data/sorce/tsi_data/daily/sorce_tsi_L3_c24h_latest.txt';
tline = 'irradiance_data.txt';
if ~exist(tline,'file') % If file is not on disk, download it.
urlwrite(url,tline);
end
fid = fopen(tline,'r');
k = 1;
while 1
tline2 = fgetl(fid);
if ~ischar(tline2),break,end;
if (length(tline2) >= 4 & tline2(1:2) == '20') %Get all the data
data(k,:) = str2num(tline2);
k = k+1;
end
end
fclose(fid);
dno = datenum(2003,2,25); %Date data began to be taken in this file
dates = [dno:dno+length(data)-1]; %Create dates vector
data(data(:,5)==0,5) = NaN; %Replace any zeroes in data with NaN.
data(:,1) = floor(data(:,1)); %Trim off decimal in date
figure(1)
plot(dates,data(:,5)) % Note that missing values (NaN) are not displayed.
hold on;
grid on;
axis tight
title('Plot of SORCE TSI data, 2012-2014','FontSize',15)
xlabel('Time (years)')
ylabel('Irradiance (W/m^2)')
datetick('x','yyyy')
xlim([734869 735964]) %Constrain x-axis to size
hold off;
2. Download videos from SDO website
I used the following Matlab code to download the videos en masse.
days = 365; % # of days we are downloading/analyzing
counter = 0;
dno = datenum(2012,8,20);
for i = dno:dno + days-1
dvo = datevec(i);
counter = counter + 1;
url = sprintf('http://sdo.gsfc.nasa.gov/assets/img/dailymov/%d/%02d/%02d/%d%02d%02d_1024_0304.mp4',...
dvo(1),dvo(2),dvo(3),dvo(1),dvo(2),dvo(3));
filename = sprintf('%d%02d%02d_1024_0304.mp4',dvo(1),dvo(2),dvo(3));
outfilename = urlwrite(url,filename)
end
3. Process videos for analysis
I used the following code to take the average of the frames in each daily video, ignoring the bad frames.
days = 365; % # of days we are analyzing
dno = datenum(2012,1,1); %Start date
meanSum = zeros(365,1); %Initialize meanSum
for t = dno:dno+days-1
dvo = datevec(t);
filename = sprintf('%d%02d%02d_1024_0304.mp4',dvo(1),dvo(2),dvo(3));
obj = VideoReader(filename);
vid = read(obj);
frames = obj.NumberOfFrames;
counter = frames;
sumImage = double(vid(:,:,:,1)); % Initialize to first image.
for k=2:frames % Read in remaining images.
rgbImage = double(vid(:,:,:,k));
mini = min(rgbImage((512-400):(512+400),(512-400):(512+400),:));
if mini == 0 %If there any black pixels in a big square in the center
counter = counter - 1; %Ignore frame
else
sumImage = sumImage + rgbImage;
end
end;
meanImage = sumImage/counter;
meanSum(ct) = sum(meanImage(:)); %Do not erase this vector; need it for plot in HW8d.m
%This text will eventually go to a log file
if counter ~= frames
fprintf('Ignored %d frames in %s\n',frames - counter,filename)
end
imwrite(uint8(meanImage),sprintf('%d%02d%02d_1024_0304.png',dvo(1),dvo(2),dvo(3)))
end
This code produces images that look like this:
I used the following Matlab code to download the videos en masse.
days = 365; % # of days we are downloading/analyzing
counter = 0;
dno = datenum(2012,8,20);
for i = dno:dno + days-1
dvo = datevec(i);
counter = counter + 1;
url = sprintf('http://sdo.gsfc.nasa.gov/assets/img/dailymov/%d/%02d/%02d/%d%02d%02d_1024_0304.mp4',...
dvo(1),dvo(2),dvo(3),dvo(1),dvo(2),dvo(3));
filename = sprintf('%d%02d%02d_1024_0304.mp4',dvo(1),dvo(2),dvo(3));
outfilename = urlwrite(url,filename)
end
3. Process videos for analysis
I used the following code to take the average of the frames in each daily video, ignoring the bad frames.
days = 365; % # of days we are analyzing
dno = datenum(2012,1,1); %Start date
meanSum = zeros(365,1); %Initialize meanSum
for t = dno:dno+days-1
dvo = datevec(t);
filename = sprintf('%d%02d%02d_1024_0304.mp4',dvo(1),dvo(2),dvo(3));
obj = VideoReader(filename);
vid = read(obj);
frames = obj.NumberOfFrames;
counter = frames;
sumImage = double(vid(:,:,:,1)); % Initialize to first image.
for k=2:frames % Read in remaining images.
rgbImage = double(vid(:,:,:,k));
mini = min(rgbImage((512-400):(512+400),(512-400):(512+400),:));
if mini == 0 %If there any black pixels in a big square in the center
counter = counter - 1; %Ignore frame
else
sumImage = sumImage + rgbImage;
end
end;
meanImage = sumImage/counter;
meanSum(ct) = sum(meanImage(:)); %Do not erase this vector; need it for plot in HW8d.m
%This text will eventually go to a log file
if counter ~= frames
fprintf('Ignored %d frames in %s\n',frames - counter,filename)
end
imwrite(uint8(meanImage),sprintf('%d%02d%02d_1024_0304.png',dvo(1),dvo(2),dvo(3)))
end
This code produces images that look like this:
The theory is that if we only consider a thin region around the edge of the sun, it will produce more reliable results. To do this, we can delete the values of the pixels in a circle using the following code.
days = 365; % # of days we are analyzing
dno = datenum(2013,1,1);
counter =0;
for t = dno:dno+days-1
counter = counter + 1;
dvo = datevec(t);
obj = imread(sprintf('%d%02d%02d_1024_0304.png',dvo(1),dvo(2),dvo(3)));
for i = 1:1024
for j = 1:1024
if sqrt((i-512)^2 + (j-512)^2) < 450
obj(i,j,:) = 0;
end
end
end
obj(200:1024,:) = 0;
A(counter) = sum(obj(:));
end
A(356)=NaN; %One missing video
B=data(3599:4329,5); %2013 and 2014's irradiance data
A_norm = normalize(A);
B_norm = normalize(B);
clf;
plot(A_norm,'r-')
hold on;
plot(B_norm,'b-')
function out=normalize(A)
I = find(isnan(A)==0);
m = mean(A(I));
s = std(A(I));
out = (A - m)/s;
days = 365; % # of days we are analyzing
dno = datenum(2013,1,1);
counter =0;
for t = dno:dno+days-1
counter = counter + 1;
dvo = datevec(t);
obj = imread(sprintf('%d%02d%02d_1024_0304.png',dvo(1),dvo(2),dvo(3)));
for i = 1:1024
for j = 1:1024
if sqrt((i-512)^2 + (j-512)^2) < 450
obj(i,j,:) = 0;
end
end
end
obj(200:1024,:) = 0;
A(counter) = sum(obj(:));
end
A(356)=NaN; %One missing video
B=data(3599:4329,5); %2013 and 2014's irradiance data
A_norm = normalize(A);
B_norm = normalize(B);
clf;
plot(A_norm,'r-')
hold on;
plot(B_norm,'b-')
function out=normalize(A)
I = find(isnan(A)==0);
m = mean(A(I));
s = std(A(I));
out = (A - m)/s;
4. Filtering and normalizing data
We can use the filter() function in Matlab to smooth the data. Below is a rough example.
z=data(3599:4329,5);
y=importdata('Sun_data.txt');
windowSize = 10;
b = (1/windowSize)*ones(1,windowSize);
a=1;
x3 = filter(b,a,z(windowSize:end));
plot(normalize(x3(windowSize:end)),'b-')
hold on
plot(normalize(y(windowSize:end)),'r-')
legend('TSI data (filtered with windowsize 10)','Luminosity data (unfiltered)')
xlabel('Days')
ylabel('Normalized intensity (unitless)')
hold off
We can use the filter() function in Matlab to smooth the data. Below is a rough example.
z=data(3599:4329,5);
y=importdata('Sun_data.txt');
windowSize = 10;
b = (1/windowSize)*ones(1,windowSize);
a=1;
x3 = filter(b,a,z(windowSize:end));
plot(normalize(x3(windowSize:end)),'b-')
hold on
plot(normalize(y(windowSize:end)),'r-')
legend('TSI data (filtered with windowsize 10)','Luminosity data (unfiltered)')
xlabel('Days')
ylabel('Normalized intensity (unitless)')
hold off