From f4b61c59dd85d6446bf23663fbe90b52147270c9 Mon Sep 17 00:00:00 2001 From: mike-dixon Date: Tue, 24 Sep 2024 17:52:04 -0600 Subject: [PATCH] Calculating spectral parameters --- .../specParams/calcSpecParams.m | 25 ++- .../specParams/plotAllMoments.m | 132 ++++++++++++++ .../specParams/plotSpecParams.m | 166 ++++++++++++++++++ .../specParams/run_specParams.m | 34 ++-- .../timeSeries/noisePeaks_smoothCorr.m | 23 ++- 5 files changed, 353 insertions(+), 27 deletions(-) create mode 100644 projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotAllMoments.m create mode 100644 projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotSpecParams.m diff --git a/projDir/qc/dataProcessing/spectralHCRproducts/specParams/calcSpecParams.m b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/calcSpecParams.m index 711d4517..c7cfd2a6 100644 --- a/projDir/qc/dataProcessing/spectralHCRproducts/specParams/calcSpecParams.m +++ b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/calcSpecParams.m @@ -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,:); @@ -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 \ No newline at end of file diff --git a/projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotAllMoments.m b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotAllMoments.m new file mode 100644 index 00000000..17e9e6fe --- /dev/null +++ b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotAllMoments.m @@ -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 \ No newline at end of file diff --git a/projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotSpecParams.m b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotSpecParams.m new file mode 100644 index 00000000..b02a2f58 --- /dev/null +++ b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/plotSpecParams.m @@ -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 \ No newline at end of file diff --git a/projDir/qc/dataProcessing/spectralHCRproducts/specParams/run_specParams.m b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/run_specParams.m index d293d7fc..d528d588 100644 --- a/projDir/qc/dataProcessing/spectralHCRproducts/specParams/run_specParams.m +++ b/projDir/qc/dataProcessing/spectralHCRproducts/specParams/run_specParams.m @@ -4,11 +4,11 @@ addpath(genpath('~/git/HCR_configuration/projDir/qc/dataProcessing/')); -project='spicule'; %socrates, aristo, cset, otrec +project='otrec'; %socrates, aristo, cset, otrec quality='ts'; %field, qc1, or qc2 -qualityCF='qc1'; +qualityCF='qc3'; freqData='10hz'; -qcVersion='v1.2'; +qcVersion='v3.2'; plotInds=0; %plotInds=(1:50:500); @@ -21,7 +21,7 @@ figdir=[dataDirCF(1:end-5),'specParams/']; -showPlot='off'; +showPlot='on'; casefile=['~/git/HCR_configuration/projDir/qc/dataProcessing/HCRproducts/caseFiles/airMotion_',project,'.txt']; @@ -301,8 +301,6 @@ for hh=1:length(dataFields1) momentsSpecSmoothCorr.(dataFields1{hh})=cat(2,momentsSpecSmoothCorr.(dataFields1{hh}),momentsSpecSmoothCorrOne.(dataFields1{hh})); end - - noiseFloorAll=cat(2,noiseFloorAll,noiseFloor); end end @@ -326,10 +324,11 @@ fieldsTs=fieldnames(momentsTime); for ll=1:length(fieldsTs) momentsTime.(fieldsTs{ll})=momentsTime.(fieldsTs{ll})(:,indTs); - momentsSpecSmoothCorr.(fieldsTs{ll})=momentsSpecSmoothCorr.(fieldsTs{ll})(:,indTs); end - - noiseFloorAll=noiseFloorAll(:,indTs); + fieldsSpec=fieldnames(momentsSpecSmoothCorr); + for ll=1:length(fieldsSpec) + momentsSpecSmoothCorr.(fieldsSpec{ll})=momentsSpecSmoothCorr.(fieldsSpec{ll})(:,indTs); + end dataCF=rmfield(dataCF,'beamWidth'); @@ -345,18 +344,9 @@ disp('Plotting ...'); momentsTime.asl=HCRrange2asl(momentsTime.range,momentsTime.elevation,momentsTime.altitude); - momentsSpecBasic.asl=HCRrange2asl(momentsSpecBasic.range,momentsSpecBasic.elevation,momentsSpecBasic.altitude); - - plotBasicMoments(momentsTime,dataCF,plotTimeAll,figdir,project,showPlot); - plotVels(momentsTime,momentsSpecBasic,momentsSpecBasicCorrRMnoise,momentsSpecSmooth,momentsSpecSmoothCorr,... - dataCF,plotTimeAll,figdir,project,showPlot); - plotWidths(momentsTime,momentsSpecBasic,momentsSpecBasicCorrRMnoise,momentsSpecSmooth,momentsSpecSmoothCorr,... - dataCF,plotTimeAll,figdir,project,showPlot); - plotSkews(momentsTime,momentsSpecBasic,momentsSpecBasicCorrRMnoise,momentsSpecSmooth,momentsSpecSmoothCorr,... - dataCF,plotTimeAll,figdir,project,showPlot); - plotKurtosis(momentsTime,momentsSpecBasic,momentsSpecBasicCorrRMnoise,momentsSpecSmooth,momentsSpecSmoothCorr,... - dataCF,plotTimeAll,figdir,project,showPlot); - - plotNoiseFloor(momentsTime,noiseFloorAll,dataCF,plotTimeAll,figdir,project,showPlot); + momentsSpecSmoothCorr.asl=HCRrange2asl(momentsSpecSmoothCorr.range,momentsSpecSmoothCorr.elevation,momentsSpecSmoothCorr.altitude); + plotAllMoments(momentsTime,dataCF,momentsSpecSmoothCorr,figdir,project,showPlot); + plotSpecParams(momentsTime,momentsSpecSmoothCorr,figdir,project,showPlot); + end diff --git a/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothCorr.m b/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothCorr.m index d8b17e80..c7a3c3b4 100644 --- a/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothCorr.m +++ b/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothCorr.m @@ -1,4 +1,4 @@ -function [powerOrig,powerOrigRMnoise,powerSmooth,powerSmoothCorr,velOut,noiseFloorAll,peakIndsAll1,peakIndsAll2]= ... +function [powerOrig,powerOrigRMnoise,powerSmooth,powerSmoothCorr,velOut,noiseFloorAllMov,peakIndsAll1,peakIndsAll2]= ... noisePeaks_smoothCorr(specDB,velIn,data,widthC,aircVel,figdir,plotTime) powerOrig=nan(size(specDB)); powerOrigRMnoise=nan(size(specDB)); @@ -132,10 +132,29 @@ end % Create new large spectrum newSpecLarge=repmat(sigWidthCorrRMnoise,1,duplicateSpec); - + newSpecLarge=cat(2,nan(1,minIndTest(ii)),newSpecLarge); newSpecLarge(end-minIndTest(ii)+1:end)=[]; + % Remove spectra pieces that are not complete + largeMask=~isnan(newSpecLarge); + statsLM=regionprops(largeMask,'Area','PixelIdxList'); + areas=[statsLM.Area]; + ua=unique(areas); + if length(ua)>1 + countRegs=nan(1,length(ua)); + for kk=1:length(ua) + countRegs(kk)=sum(areas==ua(kk)); + end + remove1=find(countRegs==1); + if ~isempty(remove1) + for ll=1:length(remove1) + regInd=find(areas==ua(remove1(ll))); + newSpecLarge(statsLM(regInd).PixelIdxList)=nan; + end + end + end + firstPowInd=find(~isnan(newSpecLarge),1,'first'); % Peaks