diff --git a/rapida/ntl/outage.py b/rapida/ntl/outage.py index 50d950e..91f3389 100644 --- a/rapida/ntl/outage.py +++ b/rapida/ntl/outage.py @@ -1,4 +1,4 @@ - +from matplotlib import pyplot as plt from datetime import datetime, timedelta import numbers import logging @@ -148,25 +148,31 @@ async def detect_outage( level = product.split('_')[0][-2:] if 'NRT' in deliverable else product[-2:] sub_dataset_name = nasa_const.SUB_DATASETS[level] daily_data = extract_bb(image_files=local_image_files, bbox=bbox,sds_name=sub_dataset_name,progress=progress) + if len(set(daily_data.ravel())) == 1: + logger.info(f'Invalid data for product {product} for {timestamp}. Skipping...') + continue positive = daily_data>=0 log_daily_data = np.zeros_like(daily_data) log_daily_data[positive] = np.log1p(daily_data[positive]) - log_daily_data = np.log1p(daily_data) daily_data_label = f'{product}_{timestamp}' arrays[daily_data_label] = log_daily_data if mask_clouds: + print('MC', product) qf_array = extract_bb(image_files=local_image_files, sds_name='QF_Cloud_Mask', bbox=bbox, progress=progress).astype('u2') cloud_confidence = (qf_array >> 6) & 0b11 is_cloudy = cloud_confidence == 3 arrays['CLOUD_MASK'] = is_cloudy - analysis_mask |= is_cloudy - arrays[f'{daily_data_label}_MASK'] = analysis_mask + product_analysis_mask = (analysis_mask | is_cloudy).astype(bool) + + else: + product_analysis_mask = analysis_mask + arrays[f'{daily_data_label}_MASK'] = product_analysis_mask log_difference, zscore, outage = utils.logdiff_outage( log_monthly_data=log_monthly_data,log_daily_data=log_daily_data, - analysis_mask=analysis_mask,percentage_drop=percentage_drop + analysis_mask=product_analysis_mask,percentage_drop=percentage_drop ) arrays[f'{daily_data_label}_LOGDIFF'] = log_difference @@ -175,9 +181,11 @@ async def detect_outage( # --- 5. UNIFIED DISPLAY & EXPORT --- #file_name = utils.get_custom_bbox_label(bbox) + file_name = get_best_semantic_label(bbox=bbox) outage_tif_path = os.path.join(dst_dir, f'{deliverable}_{file_name}.tif') outage_gpkg_path = os.path.join(dst_dir, f'{deliverable}_{file_name}.gpkg') + logger.info(f'Writing outage results to {outage_tif_path}') write_outage_tif(src_arrays=arrays, gt=gt, dst_path=outage_tif_path) if pop_vars: diff --git a/rapida/ntl/vis.py b/rapida/ntl/vis.py index 8f50978..4a64ae8 100644 --- a/rapida/ntl/vis.py +++ b/rapida/ntl/vis.py @@ -155,35 +155,51 @@ def display2(data=dict(), interpolation='nearest', title='', max_discrete_vals=5 ax.set_aspect('equal') # --- THE NEW DYNAMIC LOGIC --- - # Filter out NaNs to safely check unique values - valid_vals = a[~np.isnan(a)] + # 1. Safely filter NaNs + if np.issubdtype(a.dtype, np.floating): + valid_vals = a[~np.isnan(a)] + else: + valid_vals = a + unique_vals = np.unique(valid_vals) - # Check if this is a discrete/classification map (e.g., grid_health or masks) - if len(unique_vals) <= max_discrete_vals: - # It's discrete! Use specific colors and a custom legend. - # (Matches your Black/Yellow/Red layout for 3 values) + # 2. Handle empty (all NaNs) or uniform (all 1s) arrays natively + if len(unique_vals) < 2: + cmap = 'viridis' if 'Mask' in iname else 'magma' + # Lock the scale to 0-1 for masks so 'all 1s' doesn't default to purple + vmin, vmax = (0, 1) if ('Mask' in iname or a.dtype == bool) else (None, None) + + im = ax.imshow(a, interpolation='none', cmap=cmap, vmin=vmin, vmax=vmax) + + # Only draw a colorbar if there is at least one valid value to label + if len(unique_vals) == 1: + divider = make_axes_locatable(ax) + cax = divider.append_axes("right", size="5%", pad=0.1) + plt.colorbar(im, cax=cax, ticks=unique_vals) + + # 3. Check if this is a discrete/classification map (2 to max_discrete_vals) + elif len(unique_vals) <= max_discrete_vals: + display_a = np.zeros_like(a, dtype=int) + for idx, val in enumerate(unique_vals): + display_a[a == val] = idx + discrete_palette = ['black', 'yellow', 'red', 'cyan', 'magenta'] colors = discrete_palette[:len(unique_vals)] cmap = ListedColormap(colors) - # BoundaryNorm ensures the colorbar is split perfectly by the number of classes bounds = np.arange(len(unique_vals) + 1) - 0.5 norm = BoundaryNorm(bounds, cmap.N) - # Force interpolation='none' so classes don't blur at the edges - im = ax.imshow(a, interpolation='none', cmap=cmap, norm=norm) + im = ax.imshow(display_a, interpolation='none', cmap=cmap, norm=norm) divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.1) - # Build a discrete colorbar with centered ticks cbar = plt.colorbar(im, cax=cax, ticks=np.arange(len(unique_vals))) - # Label the ticks with their actual array values cbar.ax.set_yticklabels([f'Val: {v}' for v in unique_vals]) else: - # It's continuous! (e.g., raw radiance or log_diff) + # 4. Continuous maps cmap = 'viridis' if 'Mask' in iname else 'magma' im = ax.imshow(a, interpolation=interpolation, cmap=cmap) diff --git a/uv.lock b/uv.lock index 28321ff..2c90b2c 100644 --- a/uv.lock +++ b/uv.lock @@ -3115,6 +3115,7 @@ dependencies = [ { name = "nest-asyncio" }, { name = "netcdf4", version = "1.7.3", source = { registry = "https://pypi.org/simple" }, marker = "platform_machine == 'ARM64' and sys_platform == 'win32'" }, { name = "netcdf4", version = "1.7.4", source = { registry = "https://pypi.org/simple" }, marker = "platform_machine != 'ARM64' or sys_platform != 'win32'" }, + { name = "numpy" }, { name = "obstore" }, { name = "osm2geojson" }, { name = "overturemaps" }, @@ -3173,6 +3174,7 @@ requires-dist = [ { name = "msal" }, { name = "nest-asyncio" }, { name = "netcdf4", specifier = ">=1.7.3" }, + { name = "numpy", specifier = "<2.0.0" }, { name = "obstore", specifier = ">=0.9.5" }, { name = "osm2geojson" }, { name = "overturemaps" },