From 3ffee32c6b6bc092d9ea74200de33494ba0ed891 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 19 Dec 2023 13:59:20 -0500 Subject: [PATCH] replaced deprecated Matlab functions "nanmean", "nanstd", and "nansum" with "mean", "std", and "sum"; avoids dependency on Stats toolbox (various *.m) --- .../Create_vegopacity_8day_clim.m | 6 +-- .../get_L2_RTM_constants_tile_data.m | 4 +- .../get_model_and_obs_clim_stats.m | 48 +++++++++---------- ...get_model_and_obs_clim_stats_latlon_grid.m | 28 +++++------ .../climatology/get_model_clim_stats.m | 10 ++-- .../LDAS_App/util/postproc/write_smapL4SMqa.m | 1 - 6 files changed, 48 insertions(+), 49 deletions(-) diff --git a/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m b/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m index 3c243b27..13272941 100644 --- a/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m +++ b/src/Applications/LDAS_App/util/inputs/mwRTM_params/Create_vegopacity_8day_clim.m @@ -241,7 +241,7 @@ header = [y1 m1 d1 0 0 0 y2 m2 d2 0 0 0 tc.N_tile 1]; - tile_data = nanmean(L2_tau_tile([nidx_pre nidx nidx_nxt],:),1); + tile_data = mean(L2_tau_tile([nidx_pre nidx nidx_nxt],:),1,"omitnan"); tile_data(isnan(tile_data)) = 1.e15; fwrite( ifp, 14*4, int_precision ); % fortran_tag @@ -260,7 +260,7 @@ end % ======================== -% The final vegopacity.bin file contains data averaged (nanmean) of Asc +% The final vegopacity.bin file contains data averaged (mean) of Asc % and Des tau climatology. VOD is climatology,the other 2 parameters are % time constant with maximum spatial coverage data_clim_tile = NaN + ones(48, tc.N_tile,2); @@ -295,7 +295,7 @@ data_clim_tile(data_clim_tile < 0.) = 0.; % set small negative values to 0 % averaging A, D values -tile_data = nanmean(data_clim_tile,3); +tile_data = mean(data_clim_tile,3,"omitnan"); fname_out = strrep(fname, L2_Ascdes,'_AD_'); if fill_small_gaps diff --git a/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m b/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m index c717b9be..3125e892 100644 --- a/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m +++ b/src/Applications/LDAS_App/util/inputs/mwRTM_params/get_L2_RTM_constants_tile_data.m @@ -4,7 +4,7 @@ L2_version,start_time, end_time) % function to compute the 2 RTM constant variables: Albedo and Roughness(H) based on the -% preprocessed data. Constants are taken as the long term temporal nanmean +% preprocessed data. Constants are taken as the long term temporal mean % with maximum coverage across Asc/Desc passes. % Q. Liu 18 Jul 2022 @@ -212,7 +212,7 @@ data_clim_tile(data_clim_tile < 0.) = 0.; % set small negative values to 0 % averaging A, D values - tile_data = nanmean(data_clim_tile,2); + tile_data = mean(data_clim_tile,2,"omitnan"); eval(['out_mwRTM.',out_Para{iPara},'=tile_data;']) end diff --git a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m index 3a5cd478..9d0a7ebb 100644 --- a/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m +++ b/src/Applications/LDAS_App/util/inputs/obs_scaling_params/get_model_and_obs_clim_stats.m @@ -500,17 +500,17 @@ % and across years! % some obs can be found at multiple hours within a day % e.g. at the poles. - % **nansum of NaN's** result in zero, this need to be + % **sum(...,"omitnan") of NaNs** results in zero, this need to be % taken care of - o_data( pol(1),obs_i,angle_i,count) = nansum([o_data( pol(1),obs_i,angle_i,count); obs_obs_i' ]); - m_data( pol(1),obs_i,angle_i,count) = nansum([m_data( pol(1),obs_i,angle_i,count); obs_fcst_i']); + o_data( pol(1),obs_i,angle_i,count) = sum([o_data( pol(1),obs_i,angle_i,count); obs_obs_i' ], "omitnan"); + m_data( pol(1),obs_i,angle_i,count) = sum([m_data( pol(1),obs_i,angle_i,count); obs_fcst_i' ], "omitnan"); % X^2 - o_data2(pol(1),obs_i,angle_i,count) = nansum([o_data2(pol(1),obs_i,angle_i,count); obs_obs_i'.^2 ]); - m_data2(pol(1),obs_i,angle_i,count) = nansum([m_data2(pol(1),obs_i,angle_i,count); obs_fcst_i'.^2]); + o_data2(pol(1),obs_i,angle_i,count) = sum([o_data2(pol(1),obs_i,angle_i,count); obs_obs_i'.^2 ], "omitnan"); + m_data2(pol(1),obs_i,angle_i,count) = sum([m_data2(pol(1),obs_i,angle_i,count); obs_fcst_i'.^2 ], "omitnan"); % Sum of obs or model elements at each location - N_data(pol(1),obs_i,angle_i,count) = nansum( [N_data( pol(1),obs_i,angle_i,count); ~isnan([obs_obs_i])']); + N_data(pol(1), obs_i,angle_i,count) = sum([N_data( pol(1),obs_i,angle_i,count); ~isnan([obs_obs_i])'], "omitnan"); else @@ -523,19 +523,19 @@ % Sum of X o_data(pol(1),s_eff,angle_i,count) = ... - nansum([o_data( pol(1),s_eff,angle_i,count); repmat(obs_obs_i( i_ind), 1,length(s_eff))]); + sum([o_data( pol(1),s_eff,angle_i,count); repmat( obs_obs_i( i_ind), 1,length(s_eff))], "omitnan"); m_data(pol(1),s_eff,angle_i,count) = ... - nansum([m_data( pol(1),s_eff,angle_i,count); repmat(obs_fcst_i(i_ind), 1,length(s_eff))]); + sum([m_data( pol(1),s_eff,angle_i,count); repmat( obs_fcst_i(i_ind), 1,length(s_eff))], "omitnan"); % Sum of X^2 o_data2(pol(1),s_eff,angle_i,count) = ... - nansum([o_data2(pol(1),s_eff,angle_i,count); repmat(obs_obs_i( i_ind).^2,1,length(s_eff))]); + sum([o_data2(pol(1),s_eff,angle_i,count); repmat( obs_obs_i( i_ind).^2,1,length(s_eff))], "omitnan"); m_data2(pol(1),s_eff,angle_i,count) = ... - nansum([m_data2(pol(1),s_eff,angle_i,count); repmat(obs_fcst_i(i_ind).^2,1,length(s_eff))]); + sum([m_data2(pol(1),s_eff,angle_i,count); repmat( obs_fcst_i(i_ind).^2,1,length(s_eff))], "omitnan"); % Sum of obs or model elements at each location N_data(pol(1),s_eff,angle_i,count) = ... - nansum([N_data( pol(1),s_eff,angle_i,count); repmat(~isnan([obs_obs_i(i_ind)]),1,length(s_eff)) ]); + sum([N_data( pol(1),s_eff,angle_i,count); repmat(~isnan([obs_obs_i(i_ind)]), 1,length(s_eff))], "omitnan"); end @@ -578,41 +578,41 @@ pp = pol*Nf; - N_hscale_window = nansum(N_data(1+pol,:,:,1:w_days),4); + N_hscale_window = sum(N_data(1+pol,:,:,1:w_days), 4,"omitnan"); if w_days == 95 - N_hscale_inner_window = nansum(N_data(1+pol,:,:,((w_days+1)/2-15):((w_days+1)/2+15)),4); + N_hscale_inner_window = sum(N_data(1+pol,:,:,((w_days+1)/2-15):((w_days+1)/2+15)),4,"omitnan"); end % OBSERVATIONS %---------------- % o_data is a sum over neighbouring obs above; % here then take a sum over the time steps in the window - data_out(1+pp,:,:) = nansum(o_data(1+pol,:,:,1:w_days),4); + data_out(1+pp,:,:) = sum( o_data( 1+pol,:,:,1:w_days),4,"omitnan"); % then make the average, by dividing over the sum of the number of % timesteps and influencing obs at each location - data_out(1+pp,:,:) = data_out(1+pp,:,:)./N_hscale_window; + data_out(1+pp,:,:) = data_out( 1+pp, :,:)./N_hscale_window; %stdv_H = sqrt(E[X^2] - E[X]^2) - data_out(2+pp,:,:) = nansum(o_data2(1+pol,:,:,1:w_days),4); - data_out(2+pp,:,:) = data_out(2+pp,:,:)./N_hscale_window; - data_out(2+pp,:,:) = sqrt( data_out(2+pp,:,:) - data_out(1+pp,:,:).^2); + data_out(2+pp,:,:) = sum( o_data2( 1+pol,:,:,1:w_days),4,"omitnan"); + data_out(2+pp,:,:) = data_out( 2+pp, :,:)./N_hscale_window; + data_out(2+pp,:,:) = sqrt( data_out( 2+pp, :,:) - data_out(1+pp,:,:).^2); % MODEL %---------------- - data_out(3+pp,:,:) = nansum(m_data(1+pol,:,:,1:w_days),4); - data_out(3+pp,:,:) = data_out(3+pp,:,:)./N_hscale_window; + data_out(3+pp,:,:) = sum( m_data( 1+pol,:,:,1:w_days),4,"omitnan"); + data_out(3+pp,:,:) = data_out( 3+pp, :,:)./N_hscale_window; - data_out(4+pp,:,:) = nansum(m_data2(1+pol,:,:,1:w_days),4); - data_out(4+pp,:,:) = data_out(4+pp,:,:)./N_hscale_window; - data_out(4+pp,:,:) = sqrt( data_out(4+pp,:,:) - data_out(3+pp,:,:).^2); + data_out(4+pp,:,:) = sum( m_data2( 1+pol,:,:,1:w_days),4,"omitnan"); + data_out(4+pp,:,:) = data_out( 4+pp, :,:)./N_hscale_window; + data_out(4+pp,:,:) = sqrt( data_out( 4+pp, :,:) - data_out(3+pp,:,:).^2); data_out(5+pp,:,:) = N_hscale_window; % Toss out stats that are based on too little data - data_out([1:5]+pp,N_hscale_window