Skip to content

Commit

Permalink
Calculating spectral parameters
Browse files Browse the repository at this point in the history
  • Loading branch information
@mike-dixon
mike-dixon committed Sep 24, 2024
1 parent f55166f commit f4b61c5
Show file tree
Hide file tree
Showing 5 changed files with 353 additions and 27 deletions.
25 changes: 22 additions & 3 deletions projDir/qc/dataProcessing/spectralHCRproducts/specParams/calcSpecParams.m
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
function momentsSpecSmoothCorrOne=calcSpecParams(powerS,specVel,peaks1,peaks2,noiseFloor,ii,momentsSpecSmoothCorrOne)
loopInds=find(any(~isnan(powerS),2));

% Find the major peak if there are two
maxPeak=peaks1;
peakDiff=peaks2(:,2)-peaks1(:,2);
maxPeak(peakDiff>0,:)=peaks2(peakDiff>0,:);
Expand All @@ -11,10 +12,28 @@
powThis=powerS(jj,:);
velThis=specVel(jj,:);

rind=find(isnan(powThis),1)-1;
% Find right edge index
rind=find(~isnan(powThis),1,'last');
if isempty(rind)
rind=length(powThis);
end

% Left to right edge width
momentsSpecSmoothCorrOne.lrwidth(jj,ii)=velThis(rind)-velThis(1);
momentsSpecSmoothCorrOne.lslope(jj,ii)=(noiseFloor(jj)-maxPeak(jj,2))/(velThis(1)-velThis(maxPeak(jj,1)));
momentsSpecSmoothCorrOne.rslope(jj,ii)=(noiseFloor(jj)-maxPeak(jj,2))/(velThis(rind)-velThis(maxPeak(jj,1)));
% Left and right slope
if ~isnan(maxPeak(jj,1))
momentsSpecSmoothCorrOne.lslope(jj,ii)=(noiseFloor(jj)-maxPeak(jj,2))/(velThis(1)-velThis(maxPeak(jj,1)));
momentsSpecSmoothCorrOne.rslope(jj,ii)=(noiseFloor(jj)-maxPeak(jj,2))/(velThis(rind)-velThis(maxPeak(jj,1)));
end
% Left and right edge velocity
momentsSpecSmoothCorrOne.level(jj,ii)=velThis(1);
momentsSpecSmoothCorrOne.revel(jj,ii)=velThis(rind);
% Left and right peak velocity
if ~isnan(peaks1(jj,1))
momentsSpecSmoothCorrOne.lpvel(jj,ii)=velThis(peaks1(jj,1));
end
if ~isnan(peaks2(jj,1))
momentsSpecSmoothCorrOne.rpvel(jj,ii)=velThis(peaks2(jj,1));
end
end
end
132 changes: 132 additions & 0 deletions projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotAllMoments.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
function plotAllMoments(moments,cf,momentsSpec,figdir,project,showPlot)

momentsSpec.velRaw(:,moments.elevation>0)=-momentsSpec.velRaw(:,moments.elevation>0);
momentsSpec.skew(:,moments.elevation>0)=-momentsSpec.skew(:,moments.elevation>0);

aslGood=moments.asl(~isnan(moments.vel))./1000;
ylims=[0,max(aslGood)+0.5];

climsDbz=[-40,30];
climsLdr=[-35,-5];
climsVel=[-15,15];
climsWidth=[0,3];
climsSkew=[-3,3];
climsKurt=[-6,6];

col1=cat(1,[0,0,0],jet);
col2=cat(1,[0,0,0],velCols);

%% Figure
f1 = figure('Position',[200 500 1400 1000],'DefaultAxesFontSize',12,'visible',showPlot);

t = tiledlayout(3,2,'TileSpacing','tight','Padding','tight');

s1=nexttile(1);

cf.DBZ(isnan(cf.VEL_MASKED))=-999;

hold on
surf(moments.time,moments.asl./1000,cf.DBZ,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsDbz);
s1.Colormap=col1;
colorbar
grid on
box on
title('Time domain reflectivity (dBZ)')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s2=nexttile(2);

cf.LDR(isnan(cf.VEL_MASKED))=-999;
cf.LDR(isnan(cf.LDR))=-999;

hold on
surf(moments.time,moments.asl./1000,cf.LDR,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsLdr);
s2.Colormap=col1;
colorbar
grid on
box on
title('Time domain linear depolarization ratio (dB)')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s3=nexttile(3);

momentsSpec.velRaw(isnan(momentsSpec.velRaw))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.velRaw,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsVel);
s3.Colormap=col2;
colorbar
grid on
box on
title('Spectral domain velocity (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s4=nexttile(4);

momentsSpec.width(isnan(momentsSpec.width))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.width,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsWidth);
s4.Colormap=col1;
colorbar
grid on
box on
title('Spectral domain width (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s5=nexttile(5);

momentsSpec.skew(isnan(momentsSpec.skew))=-99;
momentsSpec.skew(isinf(momentsSpec.skew))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.skew,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsSkew);
s5.Colormap=col2;
colorbar
grid on
box on
title('Spectral domain skewness (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s6=nexttile(6);

momentsSpec.kurt(isnan(momentsSpec.kurt))=-99;
momentsSpec.kurt(isinf(momentsSpec.kurt))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.kurt,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsKurt);
s6.Colormap=col2;
colorbar
grid on
box on
title('Spectral domain kurtosis (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);


set(gcf,'PaperPositionMode','auto')
print(f1,[figdir,project,'_moments_',datestr(moments.time(1),'yyyymmdd_HHMMSS'),'_to_',datestr(moments.time(end),'yyyymmdd_HHMMSS')],'-dpng','-r0');
end
166 changes: 166 additions & 0 deletions projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotSpecParams.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,166 @@
function plotSpecParams(moments,momentsSpec,figdir,project,showPlot)

lslopeOrig=momentsSpec.lslope;
rslopeOrig=momentsSpec.rslope;
momentsSpec.lslope(:,moments.elevation>0)=-rslopeOrig(:,moments.elevation>0);
momentsSpec.rslope(:,moments.elevation>0)=-lslopeOrig(:,moments.elevation>0);

levelOrig=momentsSpec.level;
revelOrig=momentsSpec.revel;
momentsSpec.level(:,moments.elevation>0)=-revelOrig(:,moments.elevation>0);
momentsSpec.revel(:,moments.elevation>0)=-levelOrig(:,moments.elevation>0);

lpvelOrig=momentsSpec.lpvel;
rpvelOrig=momentsSpec.rpvel;
momentsSpec.lpvel(:,moments.elevation>0)=-rpvelOrig(:,moments.elevation>0);
momentsSpec.rpvel(:,moments.elevation>0)=-lpvelOrig(:,moments.elevation>0);

momentsSpec.lpvel(isnan(momentsSpec.rpvel))=nan;
momentsSpec.rpvel(isnan(momentsSpec.lpvel))=nan;

aslGood=moments.asl(~isnan(moments.vel))./1000;
ylims=[0,max(aslGood)+0.5];

climsLrwidth=[0,13];
climsLslope=[0,20];
climsRslope=[-20,0];
climsVel=[-15,15];

col1=cat(1,[0,0,0],jet);
col1r=flipud(col1);
col2=cat(1,[0,0,0],velCols);

%% Figure
f1 = figure('Position',[200 500 1400 1200],'DefaultAxesFontSize',12,'visible',showPlot);

t = tiledlayout(4,2,'TileSpacing','tight','Padding','tight');

s1=nexttile(1);

momentsSpec.lrwidth(isnan(momentsSpec.lrwidth))=-99;
momentsSpec.lrwidth(isinf(momentsSpec.lrwidth))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.lrwidth,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsLrwidth);
s1.Colormap=col1;
colorbar
grid on
box on
title('Left edge to right edge width (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s3=nexttile(3);

momentsSpec.lslope(isnan(momentsSpec.lslope))=-99;
momentsSpec.lslope(isinf(momentsSpec.lslope))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.lslope,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsLslope);
s3.Colormap=col1;
colorbar
grid on
box on
title('Left slope (dB s m^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s4=nexttile(4);

momentsSpec.rslope(isnan(momentsSpec.rslope))=99;
momentsSpec.rslope(isinf(momentsSpec.rslope))=99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.rslope,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsRslope);
s4.Colormap=col1r;
colorbar
grid on
box on
title('Right slope (dB s m^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s5=nexttile(5);

momentsSpec.level(isnan(momentsSpec.level))=-99;
momentsSpec.level(isinf(momentsSpec.level))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.level,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsVel);
s5.Colormap=col2;
colorbar
grid on
box on
title('Left edge velocity (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s6=nexttile(6);

momentsSpec.revel(isnan(momentsSpec.revel))=-99;
momentsSpec.revel(isinf(momentsSpec.revel))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.revel,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsVel);
s6.Colormap=col2;
colorbar
grid on
box on
title('Right edge velocity (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s7=nexttile(7);

momentsSpec.lpvel(isnan(momentsSpec.lpvel))=-99;
momentsSpec.lpvel(isinf(momentsSpec.lpvel))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.lpvel,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsVel);
s7.Colormap=col2;
colorbar
grid on
box on
title('Left peak velocity (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

s8=nexttile(8);

momentsSpec.rpvel(isnan(momentsSpec.rpvel))=-99;
momentsSpec.rpvel(isinf(momentsSpec.rpvel))=-99;

hold on
surf(momentsSpec.time,momentsSpec.asl./1000,momentsSpec.rpvel,'edgecolor','none');
view(2);
ylabel('Altitude (km)');
clim(climsVel);
s8.Colormap=col2;
colorbar
grid on
box on
title('Right peak velocity (m s^{-1})')
ylim(ylims);
xlim([moments.time(1),moments.time(end)]);

set(gcf,'PaperPositionMode','auto')
print(f1,[figdir,project,'_specParams_',datestr(moments.time(1),'yyyymmdd_HHMMSS'),'_to_',datestr(moments.time(end),'yyyymmdd_HHMMSS')],'-dpng','-r0');
end
Loading

0 comments on commit f4b61c5

Please sign in to comment.