import numpy as np import h5py import xarray as xr from util.util import get_grid_values_all def acspo_validate(oper_file, cspp_file, rel_tol=0.001, outfile_nc=None): h5f_oper = h5py.File(oper_file, 'r') h5f_cspp = h5py.File(cspp_file, 'r') lon_cspp = get_grid_values_all(h5f_cspp, 'lon') lat_cspp = get_grid_values_all(h5f_cspp, 'lat') print('cspp shape: ', lat_cspp.shape) c_idx = lat_cspp.shape[1] // 2 cntr_lon_cspp = lon_cspp[:, c_idx] cntr_lat_cspp = lat_cspp[:, c_idx] sst_cspp = get_grid_values_all(h5f_cspp, 'sea_surface_temperature')[0, ] l2p_flags_cspp = get_grid_values_all(h5f_cspp, 'l2p_flags')[0, ] lon_oper = get_grid_values_all(h5f_oper, 'lon') lat_oper = get_grid_values_all(h5f_oper, 'lat') cntr_lon_oper = lon_oper[:, c_idx] cntr_lat_oper = lat_oper[:, c_idx] print('oper shape: ', lat_oper.shape) sst_oper = get_grid_values_all(h5f_oper, 'sea_surface_temperature')[0, ] l2p_flags_oper = get_grid_values_all(h5f_oper, 'l2p_flags')[0, ] # generate a ndarray of boolean, True: clear or probably clear cspp_clear = (l2p_flags_cspp & (1 << 15)) == 0 oper_clear = (l2p_flags_oper & (1 << 15)) == 0 start_idx_oper, stop_idx_oper = -1, -1 start_idx_cspp, stop_idx_cspp = -1, -1 for k in range(len(cntr_lat_oper)): c_a = np.logical_and(np.isclose(cntr_lat_oper[k], cntr_lat_cspp), np.isclose(cntr_lon_oper[k], cntr_lon_cspp)) if np.size(np.nonzero(c_a)[0]) == 1: if start_idx_oper == -1: start_idx_oper = k start_idx_cspp = np.nonzero(c_a)[0][0] else: stop_idx_oper = k stop_idx_cspp = np.nonzero(c_a)[0][0] print('oper start, stop ', start_idx_oper, stop_idx_oper) print('cspp start, stop ', start_idx_cspp, stop_idx_cspp) lon_cspp = lon_cspp[start_idx_cspp:stop_idx_cspp, :] lat_cspp = lat_cspp[start_idx_cspp:stop_idx_cspp, :] lon_oper = lon_oper[start_idx_oper:stop_idx_oper, :] lat_oper = lat_oper[start_idx_oper:stop_idx_oper, :] cspp_clear = cspp_clear[start_idx_cspp:stop_idx_cspp, :] oper_clear = oper_clear[start_idx_oper:stop_idx_oper, :] overlap_shape = lon_cspp.shape print('overlap shape, size: ', overlap_shape, np.size(lon_cspp)) # check lon/lat arrays ------------------------------------ print(np.sum(np.isclose(lon_cspp, lon_oper, rtol=rel_tol))) print(np.sum(np.isclose(lat_cspp, lat_oper, rtol=rel_tol))) sst_cspp_2d = sst_cspp[start_idx_cspp:stop_idx_cspp, :] sst_oper_2d = sst_oper[start_idx_oper:stop_idx_oper, :] sst_cspp_2d = np.where(cspp_clear, sst_cspp_2d, np.nan) sst_oper_2d = np.where(oper_clear, sst_oper_2d, np.nan) sst_cspp = sst_cspp_2d.flatten() sst_oper = sst_oper_2d.flatten() oper_clear = oper_clear.flatten() cspp_clear = cspp_clear.flatten() lon_cspp = lon_cspp.flatten() lat_cspp = lat_cspp.flatten() lon_oper = lon_oper.flatten() lat_oper = lat_oper.flatten() both_clear = oper_clear & cspp_clear sst_cspp = sst_cspp[both_clear] sst_oper = sst_oper[both_clear] lon_cspp = lon_cspp[both_clear] lat_cspp = lat_cspp[both_clear] lon_oper = lon_oper[both_clear] lat_oper = lat_oper[both_clear] both_valid = np.invert(np.isnan(sst_cspp)) & np.invert(np.isnan(sst_oper)) print('number of clear and valid SSTs in both: ', np.sum(both_valid)) valid_sst_cspp = sst_cspp[both_valid] valid_sst_oper = sst_oper[both_valid] valid_lon_cspp = lon_cspp[both_valid] valid_lat_cspp = lat_cspp[both_valid] valid_lon_oper = lon_oper[both_valid] valid_lat_oper = lat_oper[both_valid] xarray_data = xr.Dataset({ 'sst_cspp': xr.DataArray(valid_sst_cspp, coords=None, dims=None, name='sst_cspp'), 'sst_oper': xr.DataArray(valid_sst_oper, coords=None, dims=None, name='sst_oper'), 'cspp_lat': xr.DataArray(valid_lat_cspp, coords=None, dims=None, name='cspp_lat'), 'cspp_lon': xr.DataArray(valid_lon_cspp, coords=None, dims=None, name='cspp_lon'), 'oper_lat': xr.DataArray(valid_lat_oper, coords=None, dims=None, name='oper_lat'), 'oper_lon': xr.DataArray(valid_lon_oper, coords=None, dims=None, name='oper_lon'), }) if outfile_nc is not None: xarray_data.to_netcdf(outfile_nc) print('fraction approx equal: ', np.sum(np.isclose(valid_sst_cspp, valid_sst_oper, rtol=rel_tol))/np.sum(both_valid)) # print(np.histogram((sst_cspp[both_valid] - sst_oper[both_valid]), bins=10)) h5f_oper.close() h5f_cspp.close() return sst_cspp_2d, sst_oper_2d def analyze_plot(): import h5py h5f_a = h5py.File('/Users/tomrink/20241107055000.nc', 'r') h5f_b = h5py.File('/Users/tomrink/20241107072000.nc', 'r') h5f_c = h5py.File('/Users/tomrink/20241107073000.nc', 'r') h5f_d = h5py.File('/Users/tomrink/20241107170000.nc', 'r') h5f_e = h5py.File('/Users/tomrink/20241107171000.nc', 'r') h5f_f = h5py.File('/Users/tomrink/20241107202000.nc', 'r') h5f_g = h5py.File('/Users/tomrink/20241107203000.nc', 'r') import numpy as np sst_cspp_all = np.concatenate([h5f_a['sst_cspp'][0, :], h5f_b['sst_cspp'][0, :], h5f_c['sst_cspp'][0, :], h5f_d['sst_cspp'][0, :], h5f_e['sst_cspp'][0, :], h5f_f['sst_cspp'][0, :], h5f_g['sst_cspp'][0, :]]) sst_cspp_all.shape sst_oper_all = np.concatenate([h5f_a['sst_oper'][0, :], h5f_b['sst_oper'][0, :], h5f_c['sst_oper'][0, :], h5f_d['sst_oper'][0, :], h5f_e['sst_oper'][0, :], h5f_f['sst_oper'][0, :], h5f_g['sst_oper'][0, :]]) sst_oper_all.shape diff = sst_cspp_all - sst_oper_all import matplotlib.pyplot as plt plt.hist(diff, bins=40) plt.yscale('log') plt.title('ACSPO SST (CSPP - OPER)') import h5py h5f_a = h5py.File('/Users/tomrink/20241107055000.nc', 'r') h5f_b = h5py.File('/Users/tomrink/20241107072000.nc', 'r') h5f_c = h5py.File('/Users/tomrink/20241107073000.nc', 'r') h5f_d = h5py.File('/Users/tomrink/20241107170000.nc', 'r') h5f_e = h5py.File('/Users/tomrink/20241107171000.nc', 'r') h5f_f = h5py.File('/Users/tomrink/20241107202000.nc', 'r') h5f_g = h5py.File('/Users/tomrink/20241107203000.nc', 'r') import numpy as np sst_cspp_all = np.concatenate([h5f_a['sst_cspp'][0, :], h5f_b['sst_cspp'][0, :], h5f_c['sst_cspp'][0, :], h5f_d['sst_cspp'][0, :], h5f_e['sst_cspp'][0, :], h5f_f['sst_cspp'][0, :], h5f_g['sst_cspp'][0, :]]) sst_oper_all = np.concatenate([h5f_a['sst_oper'][0, :], h5f_b['sst_oper'][0, :], h5f_c['sst_oper'][0, :], h5f_d['sst_oper'][0, :], h5f_e['sst_oper'][0, :], h5f_f['sst_oper'][0, :], h5f_g['sst_oper'][0, :]]) diff = sst_cspp_all - sst_oper_all import matplotlib.pyplot as plt plt.hist(diff, bins=40) plt.yscale('log') plt.xlabel('SST (\u00b0K)') plt.title('ACSPO SST (CSPP - OPER), Clear Sky, 2024-11_07') import numpy as np c_a = np.isclose(sst_cspp_all, sst_oper_all, rtol=0.001) c_a.shape sst_cspp_all.shape c_a_i = np.invert(c_a) plt.scatter(sst_cspp_all[c_a_i], diff[c_a_i]) plt.close() plt.scatter(sst_cspp_all[c_a_i], diff[c_a_i], s=1) plt.title('ACSPO SST (CSPP-OPER) > rel_tol=0.001 vs SST') plt.xlabel('SST (\u00b0K)') plt.ylabel('SST (\u00b0K)') plt.show() import h5py h5f_a = h5py.File('/Users/tomrink/20241107055000.nc', 'r') h5f_b = h5py.File('/Users/tomrink/20241107072000.nc', 'r') h5f_c = h5py.File('/Users/tomrink/20241107073000.nc', 'r') h5f_d = h5py.File('/Users/tomrink/20241107170000.nc', 'r') h5f_e = h5py.File('/Users/tomrink/20241107171000.nc', 'r') h5f_f = h5py.File('/Users/tomrink/20241107202000.nc', 'r') h5f_g = h5py.File('/Users/tomrink/20241107203000.nc', 'r') import numpy as np sst_cspp_all = np.concatenate( [h5f_a['sst_cspp'][0, :], h5f_b['sst_cspp'][0, :], h5f_c['sst_cspp'][0, :], h5f_d['sst_cspp'][0, :], h5f_e['sst_cspp'][0, :], h5f_f['sst_cspp'][0, :], h5f_g['sst_cspp'][0, :]]) sst_cspp_all.shape sst_oper_all = np.concatenate( [h5f_a['sst_oper'][0, :], h5f_b['sst_oper'][0, :], h5f_c['sst_oper'][0, :], h5f_d['sst_oper'][0, :], h5f_e['sst_oper'][0, :], h5f_f['sst_oper'][0, :], h5f_g['sst_oper'][0, :]]) sst_oper_all.shape diff = sst_cspp_all - sst_oper_all import matplotlib.pyplot as plt plt.hist(diff, bins=40) plt.yscale('log') plt.title('ACSPO SST (CSPP - OPER)') import h5py h5f_a = h5py.File('/Users/tomrink/20241107055000.nc', 'r') h5f_b = h5py.File('/Users/tomrink/20241107072000.nc', 'r') h5f_c = h5py.File('/Users/tomrink/20241107073000.nc', 'r') h5f_d = h5py.File('/Users/tomrink/20241107170000.nc', 'r') h5f_e = h5py.File('/Users/tomrink/20241107171000.nc', 'r') h5f_f = h5py.File('/Users/tomrink/20241107202000.nc', 'r') h5f_g = h5py.File('/Users/tomrink/20241107203000.nc', 'r') import numpy as np sst_cspp_all = np.concatenate( [h5f_a['sst_cspp'][0, :], h5f_b['sst_cspp'][0, :], h5f_c['sst_cspp'][0, :], h5f_d['sst_cspp'][0, :], h5f_e['sst_cspp'][0, :], h5f_f['sst_cspp'][0, :], h5f_g['sst_cspp'][0, :]]) sst_oper_all = np.concatenate( [h5f_a['sst_oper'][0, :], h5f_b['sst_oper'][0, :], h5f_c['sst_oper'][0, :], h5f_d['sst_oper'][0, :], h5f_e['sst_oper'][0, :], h5f_f['sst_oper'][0, :], h5f_g['sst_oper'][0, :]]) diff = sst_cspp_all - sst_oper_all import matplotlib.pyplot as plt plt.hist(diff, bins=40) plt.yscale('log') plt.xlabel('SST (\u00b0K)') plt.title('ACSPO SST (CSPP - OPER), Clear Sky, 2024-11_07') import numpy as np c_a = np.isclose(sst_cspp_all, sst_oper_all, rtol=0.001) c_a.shape sst_cspp_all.shape c_a_i = np.invert(c_a) plt.scatter(sst_cspp_all[c_a_i], diff[c_a_i]) plt.close() plt.scatter(sst_cspp_all[c_a_i], diff[c_a_i], s=1) plt.title('ACSPO SST (CSPP-OPER) > rel_tol=0.001 vs SST') plt.xlabel('SST (\u00b0K)') plt.ylabel('SST (\u00b0K)') plt.show()