# -*- coding: utf-8 -*- """ Sistema de Predicción de Peligro de Incendios Forestales de México Índice de Verdor Relativo Archivo de funciones adicionales """ import numpy as np import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling from rasterio.transform import from_origin import tempfile import os from typing import Optional, Union, Tuple class RasterioProcessor: """ Clase para realizar operaciones de geoprocesamiento con objetos rasterio.io.DatasetReader """ def __init__(self): self.temp_files = [] def __del__(self): """Limpia los archivos temporales al destruir la instancia""" self.cleanup() def cleanup(self): """Elimina todos los archivos temporales creados""" for temp_file in self.temp_files: try: if os.path.exists(temp_file): os.unlink(temp_file) except: pass self.temp_files = [] def _create_temp_raster(self, data: np.ndarray, profile: dict) -> rasterio.io.DatasetReader: """ Crea un raster temporal con los datos y perfil proporcionados Args: data: Array numpy con los datos del raster profile: Perfil (metadatos) del raster Returns: rasterio.io.DatasetReader: Dataset del raster temporal """ # Crear archivo temporal temp_file = tempfile.mktemp(suffix='.tif') self.temp_files.append(temp_file) return self._write_raster_to_file(data, profile, temp_file) def _write_raster_to_file(self, data: np.ndarray, profile: dict, output_path: str) -> rasterio.io.DatasetReader: """ Escribe un raster a un archivo Args: data: Array numpy con los datos del raster profile: Perfil (metadatos) del raster output_path: Ruta del archivo de salida Returns: rasterio.io.DatasetReader: Dataset del raster escrito """ # Actualizar el perfil con las dimensiones y tipo de datos profile.update({ 'height': data.shape[1] if len(data.shape) == 3 else data.shape[0], 'width': data.shape[2] if len(data.shape) == 3 else data.shape[1], 'count': 1 if len(data.shape) == 2 else data.shape[0], 'dtype': data.dtype, 'driver': 'GTiff' }) # Escribir el raster with rasterio.open(output_path, 'w', **profile) as dst: if len(data.shape) == 2: dst.write(data, 1) else: for i in range(data.shape[0]): dst.write(data[i], i + 1) # Abrir y retornar el dataset return rasterio.open(output_path) def save_to_tif(self, dataset: rasterio.io.DatasetReader, output_path: str, compress: Optional[str] = None, nodata: Optional[Union[int, float]] = None) -> rasterio.io.DatasetReader: """ Guarda un dataset en un archivo .tif Args: dataset: Dataset a guardar (rasterio.io.DatasetReader) output_path: Ruta del archivo de salida .tif compress: Opción de compresión (ej: 'LZW', 'DEFLATE', 'PACKBITS') nodata: Valor nodata para el archivo de salida Returns: rasterio.io.DatasetReader: Dataset del archivo guardado """ # Leer datos data = dataset.read() # Crear perfil de salida profile = dataset.profile.copy() profile.update({ 'driver': 'GTiff', 'count': dataset.count, 'height': dataset.height, 'width': dataset.width, 'dtype': data.dtype }) # Aplicar compresión si se especifica if compress: profile.update({ 'compress': compress, 'tiled': True, 'blockxsize': 256, 'blockysize': 256 }) # Actualizar valor nodata si se especifica if nodata is not None: profile['nodata'] = nodata return self._write_raster_to_file(data, profile, output_path) def reproject_to_match(self, src_dataset: rasterio.io.DatasetReader, target_dataset: rasterio.io.DatasetReader, output_path: Optional[str] = None, resampling_method: Resampling = Resampling.nearest) -> rasterio.io.DatasetReader: """ Reprojecta un dataset para que coincida con el extent y tamaño de pixel de otro Args: src_dataset: Dataset de origen (rasterio.io.DatasetReader) target_dataset: Dataset objetivo para copiar extent y tamaño de pixel output_path: Ruta opcional para guardar el resultado como .tif resampling_method: Método de remuestreo (por defecto: nearest) Returns: rasterio.io.DatasetReader: Dataset reproyectado """ # Obtener transformación del dataset objetivo target_crs = target_dataset.crs target_bounds = target_dataset.bounds target_width = target_dataset.width target_height = target_dataset.height # Calcular transformación por defecto para reproyección transform, width, height = calculate_default_transform( src_dataset.crs, target_crs, src_dataset.width, src_dataset.height, *src_dataset.bounds ) # Crear array para datos reproyectados data = src_dataset.read() reprojected_data = np.zeros((src_dataset.count, target_height, target_width), dtype=data.dtype) # Reproyectar reproject( source=data, destination=reprojected_data, src_transform=src_dataset.transform, src_crs=src_dataset.crs, dst_transform=target_dataset.transform, dst_crs=target_crs, resampling=resampling_method ) # Crear perfil para el raster reproyectado profile = target_dataset.profile.copy() profile.update({ 'count': src_dataset.count, 'dtype': data.dtype, 'nodata': src_dataset.nodata }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(reprojected_data, profile, output_path) else: return self._create_temp_raster(reprojected_data, profile) def multiply_datasets(self, dataset1: rasterio.io.DatasetReader, dataset2: rasterio.io.DatasetReader, output_path: Optional[str] = None, ignore_nodata: bool = True) -> rasterio.io.DatasetReader: """ Multiplica los datos de dos datasets con opción de ignorar nodata Args: dataset1: Primer dataset (rasterio.io.DatasetReader) dataset2: Segundo dataset (rasterio.io.DatasetReader) output_path: Ruta opcional para guardar el resultado como .tif ignore_nodata: Si True, ignora valores nodata en la operación Returns: rasterio.io.DatasetReader: Dataset con el resultado de la multiplicación """ # Leer datos de ambos datasets data1 = dataset1.read() data2 = dataset2.read() # Asegurar que tengan la misma forma if data1.shape != data2.shape: # Reprojectar dataset2 para que coincida con dataset1 dataset2_reproj = self.reproject_to_match(dataset2, dataset1) data2 = dataset2_reproj.read() dataset2_reproj.close() # Crear copias para modificar si es necesario data1_proc = data1.copy() data2_proc = data2.copy() # Manejar nodata si se ignora if ignore_nodata: # Manejar nodata de dataset1 if dataset1.nodata is not None: if np.issubdtype(data1.dtype, np.floating): mask1 = np.isnan(data1) if np.isnan(dataset1.nodata) else np.isclose(data1, dataset1.nodata) else: mask1 = (data1 == dataset1.nodata) # Reemplazar nodata con 1 (elemento neutro para multiplicación) data1_proc[mask1] = 1 # Manejar nodata de dataset2 if dataset2.nodata is not None: if np.issubdtype(data2.dtype, np.floating): mask2 = np.isnan(data2) if np.isnan(dataset2.nodata) else np.isclose(data2, dataset2.nodata) else: mask2 = (data2 == dataset2.nodata) # Reemplazar nodata con 1 data2_proc[mask2] = 1 # La máscara de nodata final es la combinación if dataset1.nodata is not None and dataset2.nodata is not None: final_nodata_mask = np.logical_or(mask1, mask2) elif dataset1.nodata is not None: final_nodata_mask = mask1 elif dataset2.nodata is not None: final_nodata_mask = mask2 else: final_nodata_mask = None else: final_nodata_mask = None # Realizar la multiplicación result = data1_proc * data2_proc # Aplicar nodata al resultado si es necesario if final_nodata_mask is not None and np.any(final_nodata_mask): # Determinar valor nodata result_nodata = dataset1.nodata if dataset1.nodata is not None else dataset2.nodata if result_nodata is not None: # Asignar nodata if np.isnan(result_nodata): result = result.astype(np.float32) result[final_nodata_mask] = np.nan else: if np.issubdtype(result.dtype, np.floating): result[final_nodata_mask] = float(result_nodata) else: result[final_nodata_mask] = int(result_nodata) # Determinar dtype apropiado if np.issubdtype(data1.dtype, np.floating) or np.issubdtype(data2.dtype, np.floating): result_dtype = np.float32 result = result.astype(result_dtype) else: result_dtype = result.dtype # Crear perfil de salida profile = dataset1.profile.copy() profile.update({ 'count': result.shape[0] if len(result.shape) == 3 else 1, 'dtype': result_dtype, 'nodata': dataset1.nodata if not ignore_nodata else (dataset1.nodata if dataset1.nodata is not None else dataset2.nodata) }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def multiply_by_scalar(self, dataset: rasterio.io.DatasetReader, scalar: Union[int, float], output_path: Optional[str] = None, preserve_nodata: bool = True, handle_nan: bool = True) -> rasterio.io.DatasetReader: """ Multiplica los datos de un dataset por un escalar, preservando valores nodata Args: dataset: Dataset a multiplicar (rasterio.io.DatasetReader) scalar: Número escalar para multiplicar output_path: Ruta opcional para guardar el resultado como .tif preserve_nodata: Si True, preserva los valores nodata sin multiplicar handle_nan: Si True, maneja NaN values como nodata adicionalmente Returns: rasterio.io.DatasetReader: Dataset con el resultado de la multiplicación """ # Leer datos data = dataset.read() original_dtype = data.dtype # Inicializar máscara de valores a preservar preserve_mask = np.zeros_like(data, dtype=bool) # Identificar valores a preservar if preserve_nodata: # 1. Valores nodata definidos if dataset.nodata is not None: # Manejar comparación segura para float if np.issubdtype(data.dtype, np.floating): # Para floats, usar comparación con tolerancia preserve_mask = np.isclose(data, dataset.nodata) else: # Para enteros, comparación directa preserve_mask = (data == dataset.nodata) # 2. Valores NaN (si se solicitó) if handle_nan and np.issubdtype(data.dtype, np.floating): nan_mask = np.isnan(data) preserve_mask = np.logical_or(preserve_mask, nan_mask) # Multiplicar datos result = data * scalar # Preservar valores originales donde corresponda if preserve_nodata and np.any(preserve_mask): result[preserve_mask] = data[preserve_mask] # Manejar conversión de tipos si es necesario if isinstance(scalar, float) and np.issubdtype(original_dtype, np.integer): # Escalar flotante con datos enteros -> convertir a float result = result.astype(np.float32) # Si había nodata definido, mantenerlo if preserve_nodata and dataset.nodata is not None: # Buscar posiciones originales con nodata if np.issubdtype(original_dtype, np.floating): nodata_positions = np.isclose(data, dataset.nodata) else: nodata_positions = (data == dataset.nodata) # Asignar nodata convertido a float result[nodata_positions] = float(dataset.nodata) # Crear perfil de salida profile = dataset.profile.copy() # Actualizar nodata si se convirtió a float new_nodata = dataset.nodata if isinstance(scalar, float) and np.issubdtype(original_dtype, np.integer): if dataset.nodata is not None: new_nodata = float(dataset.nodata) profile.update({ 'dtype': result.dtype, 'nodata': new_nodata }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def subtract_dataset_to_scalar(self, dataset: rasterio.io.DatasetReader, scalar: Union[int, float], output_path: Optional[str] = None) -> rasterio.io.DatasetReader: """ Resta un número escalar de un dataset Args: dataset: Dataset de entrada (rasterio.io.DatasetReader) scalar: Número escalar a restar output_path: Ruta opcional para guardar el resultado como .tif Returns: rasterio.io.DatasetReader: Dataset con el resultado de la resta """ # Leer datos y restar escalar data = dataset.read() result = data - scalar # Crear perfil de salida profile = dataset.profile.copy() profile.update({ 'dtype': result.dtype }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def subtract_datasets(self, dataset1: rasterio.io.DatasetReader, dataset2: rasterio.io.DatasetReader, output_path: Optional[str] = None, ignore_nodata: bool = False, nodata_output: Optional[Union[int, float]] = None, fill_method: str = 'propagate') -> rasterio.io.DatasetReader: """ Resta un dataset de otro (dataset1 - dataset2) Args: dataset1: Dataset del que se resta (rasterio.io.DatasetReader) dataset2: Dataset a restar (rasterio.io.DatasetReader) output_path: Ruta opcional para guardar el resultado como .tif ignore_nodata: Si True, ignora valores nodata en la operación nodata_output: Valor nodata para el resultado (si None, usa dataset1.nodata) fill_method: Método para llenar nodata en el resultado: 'propagate' - propaga nodata de cualquier dataset 'dataset1' - usa valores de dataset1 donde dataset2 tiene nodata 'dataset2' - usa valores de dataset2 donde dataset1 tiene nodata 'value' - usa un valor específico (requiere nodata_output) Returns: rasterio.io.DatasetReader: Dataset con el resultado de la resta """ # Leer datos de ambos datasets data1 = dataset1.read() data2 = dataset2.read() # Asegurar que tengan la misma forma if data1.shape != data2.shape: # Reprojectar dataset2 para que coincida con dataset1 dataset2_reproj = self.reproject_to_match(dataset2, dataset1) data2 = dataset2_reproj.read() dataset2_reproj.close() # Crear máscaras de nodata si es necesario nodata_mask1 = None nodata_mask2 = None nodata_mask_result = None if ignore_nodata: # Crear copias de los datos para modificar data1_processed = data1.copy() data2_processed = data2.copy() # Identificar valores nodata en cada dataset if dataset1.nodata is not None: if np.issubdtype(data1.dtype, np.floating): nodata_mask1 = np.isnan(data1) if np.isnan(dataset1.nodata) else np.isclose(data1, dataset1.nodata) else: nodata_mask1 = (data1 == dataset1.nodata) if dataset2.nodata is not None: if np.issubdtype(data2.dtype, np.floating): nodata_mask2 = np.isnan(data2) if np.isnan(dataset2.nodata) else np.isclose(data2, dataset2.nodata) else: nodata_mask2 = (data2 == dataset2.nodata) # Manejar nodata según el método seleccionado if fill_method == 'propagate': # Crear máscara combinada de nodata if nodata_mask1 is not None or nodata_mask2 is not None: if nodata_mask1 is None: nodata_mask_result = nodata_mask2 elif nodata_mask2 is None: nodata_mask_result = nodata_mask1 else: nodata_mask_result = np.logical_or(nodata_mask1, nodata_mask2) # Para la operación, reemplazar nodata con 0 temporalmente if nodata_mask1 is not None: data1_processed[nodata_mask1] = 0 if nodata_mask2 is not None: data2_processed[nodata_mask2] = 0 elif fill_method == 'dataset1': # Donde dataset2 tiene nodata, usar solo dataset1 (resultado = dataset1) if nodata_mask2 is not None: # Para la operación, hacer dataset2 = 0 donde tiene nodata data2_processed[nodata_mask2] = 0 # El resultado mantendrá el valor de dataset1 en esas posiciones nodata_mask_result = nodata_mask1 # Solo nodata de dataset1 elif fill_method == 'dataset2': # Donde dataset1 tiene nodata, usar solo dataset2 (resultado = -dataset2) if nodata_mask1 is not None: # Para la operación, hacer dataset1 = 0 donde tiene nodata data1_processed[nodata_mask1] = 0 # El resultado será negativo de dataset2 en esas posiciones nodata_mask_result = nodata_mask2 # Solo nodata de dataset2 elif fill_method == 'value': # Usar un valor específico para nodata en la operación fill_value = nodata_output if nodata_output is not None else 0 if nodata_mask1 is not None: data1_processed[nodata_mask1] = fill_value if nodata_mask2 is not None: data2_processed[nodata_mask2] = fill_value # La máscara de resultado sigue siendo la combinada if nodata_mask1 is not None or nodata_mask2 is not None: if nodata_mask1 is None: nodata_mask_result = nodata_mask2 elif nodata_mask2 is None: nodata_mask_result = nodata_mask1 else: nodata_mask_result = np.logical_or(nodata_mask1, nodata_mask2) else: raise ValueError(f"Método fill_method no válido: {fill_method}") # Realizar la resta con datos procesados result = data1_processed - data2_processed # Aplicar máscara de nodata al resultado if nodata_mask_result is not None: # Determinar valor nodata para el resultado if nodata_output is not None: result_nodata = nodata_output elif dataset1.nodata is not None: result_nodata = dataset1.nodata elif dataset2.nodata is not None: result_nodata = dataset2.nodata else: # Valor por defecto según el tipo de datos if np.issubdtype(result.dtype, np.floating): result_nodata = np.nan elif np.issubdtype(result.dtype, np.integer): result_nodata = -9999 else: result_nodata = None # Asignar nodata al resultado if result_nodata is not None: # Convertir si es necesario if np.isnan(result_nodata): result = result.astype(np.float32) result[nodata_mask_result] = np.nan else: # Asegurar tipo de datos compatible if np.issubdtype(result.dtype, np.floating): result[nodata_mask_result] = float(result_nodata) else: result[nodata_mask_result] = int(result_nodata) else: # Comportamiento original: incluir nodata en la operación result = data1 - data2 nodata_mask_result = None # Determinar dtype apropiado para el resultado # Si hay mezcla de tipos, usar el más preciso if data1.dtype != data2.dtype: if np.issubdtype(data1.dtype, np.floating) or np.issubdtype(data2.dtype, np.floating): result_dtype = np.float32 result = result.astype(result_dtype) else: result_dtype = result.dtype else: result_dtype = result.dtype # Crear perfil de salida profile = dataset1.profile.copy() # Determinar valor nodata para el perfil if ignore_nodata: if nodata_output is not None: profile_nodata = nodata_output elif dataset1.nodata is not None: profile_nodata = dataset1.nodata else: profile_nodata = None else: profile_nodata = dataset1.nodata profile.update({ 'count': result.shape[0] if len(result.shape) == 3 else 1, 'dtype': result_dtype, 'nodata': profile_nodata }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def add_scalar(self, dataset: rasterio.io.DatasetReader, scalar: Union[int, float], output_path: Optional[str] = None) -> rasterio.io.DatasetReader: """ Suma un número escalar a un dataset Args: dataset: Dataset de entrada (rasterio.io.DatasetReader) scalar: Número escalar a sumar output_path: Ruta opcional para guardar el resultado como .tif Returns: rasterio.io.DatasetReader: Dataset con el resultado de la suma """ # Leer datos y sumar escalar data = dataset.read() result = data + scalar # Crear perfil de salida profile = dataset.profile.copy() profile.update({ 'dtype': result.dtype }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def add_datasets(self, dataset1: rasterio.io.DatasetReader, dataset2: rasterio.io.DatasetReader, output_path: Optional[str] = None) -> rasterio.io.DatasetReader: """ Suma dos datasets Args: dataset1: Primer dataset (rasterio.io.DatasetReader) dataset2: Segundo dataset (rasterio.io.DatasetReader) output_path: Ruta opcional para guardar el resultado como .tif Returns: rasterio.io.DatasetReader: Dataset con el resultado de la suma """ # Leer datos de ambos datasets data1 = dataset1.read() data2 = dataset2.read() # Asegurar que tengan la misma forma if data1.shape != data2.shape: # Reprojectar dataset2 para que coincida con dataset1 dataset2_reproj = self.reproject_to_match(dataset2, dataset1) data2 = dataset2_reproj.read() dataset2_reproj.close() # Sumar los datos result = data1 + data2 # Crear perfil de salida profile = dataset1.profile.copy() profile.update({ 'count': result.shape[0] if len(result.shape) == 3 else 1, 'dtype': result.dtype }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def set_nodata_value(self, dataset: rasterio.io.DatasetReader, nodata_value: Union[int, float], output_path: Optional[str] = None) -> rasterio.io.DatasetReader: """ Asigna un nuevo valor para los píxeles nodata Args: dataset: Dataset de entrada (rasterio.io.DatasetReader) nodata_value: Nuevo valor para nodata output_path: Ruta opcional para guardar el resultado como .tif Returns: rasterio.io.DatasetReader: Dataset con el cambio aplicado en nodata """ # Leer datos data = dataset.read() # Crear perfil con el nuevo valor nodata profile = dataset.profile.copy() profile.update({ 'nodata': nodata_value }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(data, profile, output_path) else: return self._create_temp_raster(data, profile) def float32_to_int16(self, dataset: rasterio.io.DatasetReader, output_path: Optional[str] = None, nodata_value: int = -32768, clamp_values: bool = True) -> rasterio.io.DatasetReader: """ Conversión simple de float32 a int16 con manejo básico. Args: dataset: Dataset de entrada en float32 output_path: Ruta opcional para guardar el resultado nodata_value: Valor nodata para el output clamp_values: Si True, recorta valores fuera del rango int16 Returns: rasterio.io.DatasetReader: Dataset convertido a int16 """ # Verificar tipo de datos if dataset.dtypes[0] != 'float32': raise ValueError(f"Se requiere float32, pero el dataset es {dataset.dtypes[0]}") # Leer datos data = dataset.read() # Manejar nodata original original_nodata = dataset.nodata nodata_mask = None if original_nodata is not None: if np.isnan(original_nodata): nodata_mask = np.isnan(data) else: nodata_mask = np.isclose(data, original_nodata) # Convertir directamente (puede perder precisión) result = data.astype(np.int16) # Recortar valores si se solicita if clamp_values: result = np.clip(result, -32768, 32767) # Manejar nodata if nodata_mask is not None and np.any(nodata_mask): result[nodata_mask] = nodata_value # Crear perfil de salida profile = dataset.profile.copy() profile.update({ 'dtype': 'int16', 'nodata': nodata_value }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def divide_datasets(self, dataset1: rasterio.io.DatasetReader, dataset2: rasterio.io.DatasetReader, output_path: Optional[str] = None, safe_division: bool = True) -> rasterio.io.DatasetReader: """ División simple de dos datasets con manejo seguro por defecto Args: dataset1: Dataset numerador dataset2: Dataset denominador output_path: Ruta opcional para guardar el resultado safe_division: Si True, evita división por cero y maneja nodata Returns: rasterio.io.DatasetReader: Dataset con el resultado de la división """ # Leer datos data1 = dataset1.read() data2 = dataset2.read() # Asegurar misma forma if data1.shape != data2.shape: dataset2_reproj = self.reproject_to_match(dataset2, dataset1) data2 = dataset2_reproj.read() dataset2_reproj.close() # Convertir a float para división data1_float = data1.astype(np.float32) data2_float = data2.astype(np.float32) if safe_division: # Evitar división por cero data2_float[data2_float == 0] = np.nan # Manejar nodata if dataset1.nodata is not None: if np.issubdtype(data1.dtype, np.floating): nodata_mask1 = np.isnan(data1) if np.isnan(dataset1.nodata) else np.isclose(data1, dataset1.nodata) else: nodata_mask1 = (data1 == dataset1.nodata) data1_float[nodata_mask1] = np.nan if dataset2.nodata is not None: if np.issubdtype(data2.dtype, np.floating): nodata_mask2 = np.isnan(data2) if np.isnan(dataset2.nodata) else np.isclose(data2, dataset2.nodata) else: nodata_mask2 = (data2 == dataset2.nodata) data2_float[nodata_mask2] = np.nan # Realizar división with np.errstate(divide='ignore', invalid='ignore'): result = np.divide(data1_float, data2_float) # Reemplazar NaN con nodata result_nan_mask = np.isnan(result) if np.any(result_nan_mask): result[result_nan_mask] = np.nan # Crear perfil profile = dataset1.profile.copy() profile.update({ 'count': result.shape[0] if len(result.shape) == 3 else 1, 'dtype': 'float32', 'nodata': np.nan }) # Guardar si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def divide_by_scalar(self, dataset: rasterio.io.DatasetReader, scalar: Union[int, float], output_path: Optional[str] = None, ignore_nodata: bool = True) -> rasterio.io.DatasetReader: """ Divide los datos de un dataset por un escalar Args: dataset: Dataset a dividir (rasterio.io.DatasetReader) scalar: Número escalar para dividir output_path: Ruta opcional para guardar el resultado como .tif ignore_nodata: Si True, ignora valores nodata en la operación Returns: rasterio.io.DatasetReader: Dataset con el resultado de la división """ # Verificar que el escalar no sea 0 if scalar == 0: raise ValueError("No se puede dividir por cero. Use un valor diferente a 0.") # Leer datos del dataset data = dataset.read() # Crear copia para modificar si es necesario data_proc = data.copy() # Manejar nodata si se ignora if ignore_nodata and dataset.nodata is not None: # Crear máscara para valores nodata if np.issubdtype(data.dtype, np.floating): nodata_mask = np.isnan(data) if np.isnan(dataset.nodata) else np.isclose(data, dataset.nodata) else: nodata_mask = (data == dataset.nodata) # Para la operación, reemplazar nodata con 0 temporalmente # (0 dividido por cualquier cosa sigue siendo 0) data_proc[nodata_mask] = 0 else: nodata_mask = None # Realizar la división result = data_proc / scalar # Restaurar valores nodata si se ignoraron if ignore_nodata and nodata_mask is not None and np.any(nodata_mask): # Usar el valor nodata original result_nodata = dataset.nodata # Asignar nodata al resultado if result_nodata is not None: if np.isnan(result_nodata): # Para NaN, asegurar que el resultado sea float result = result.astype(np.float32) result[nodata_mask] = np.nan else: # Para otros valores, convertir si es necesario if np.issubdtype(result.dtype, np.floating): result[nodata_mask] = float(result_nodata) else: result[nodata_mask] = int(result_nodata) # Determinar dtype apropiado para el resultado # Si el escalar es entero y los datos son enteros, mantener int si es posible # De lo contrario, convertir a float if isinstance(scalar, int) and np.issubdtype(data.dtype, np.integer): # Verificar si la división produce valores enteros # Para simplificar, convertimos a float de todos modos para mayor precisión result_dtype = np.float32 result = result.astype(result_dtype) else: # Para división con float o datos float, usar float32 result_dtype = np.float32 result = result.astype(result_dtype) # Crear perfil de salida profile = dataset.profile.copy() # Determinar valor nodata para el perfil if ignore_nodata: profile_nodata = dataset.nodata else: profile_nodata = dataset.nodata profile.update({ 'dtype': result_dtype, 'nodata': profile_nodata }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile) def subtract_scalar_from_dataset(self, dataset: rasterio.io.DatasetReader, scalar: Union[int, float], output_path: Optional[str] = None, ignore_nodata: bool = True) -> rasterio.io.DatasetReader: """ Resta un número escalar de un dataset Args: dataset: Dataset de entrada (rasterio.io.DatasetReader) scalar: Número escalar a restar output_path: Ruta opcional para guardar el resultado como .tif ignore_nodata: Si True, ignora valores nodata en la operación Returns: rasterio.io.DatasetReader: Dataset con el resultado de la resta """ # Leer datos del dataset data = dataset.read() # Crear copia para modificar si es necesario data_proc = data.copy() # Manejar nodata si se ignora if ignore_nodata and dataset.nodata is not None: # Crear máscara para valores nodata if np.issubdtype(data.dtype, np.floating): nodata_mask = np.isnan(data) if np.isnan(dataset.nodata) else np.isclose(data, dataset.nodata) else: nodata_mask = (data == dataset.nodata) # Para la operación, reemplazar nodata con 0 temporalmente # (0 - escalar = -escalar, pero luego restaurará el nodata) data_proc[nodata_mask] = 0 else: nodata_mask = None # Realizar la resta #result = data_proc - scalar result = scalar - data_proc # Restaurar valores nodata si se ignoraron if ignore_nodata and nodata_mask is not None and np.any(nodata_mask): # Usar el valor nodata original result_nodata = dataset.nodata # Asignar nodata al resultado if result_nodata is not None: if np.isnan(result_nodata): # Para NaN, asegurar que el resultado sea float result = result.astype(np.float32) result[nodata_mask] = np.nan else: # Para otros valores, convertir si es necesario if np.issubdtype(result.dtype, np.floating): result[nodata_mask] = float(result_nodata) else: result[nodata_mask] = int(result_nodata) # Determinar dtype apropiado para el resultado # Si el escalar es float y los datos son enteros, convertir a float if isinstance(scalar, float) and np.issubdtype(data.dtype, np.integer): result_dtype = np.float32 result = result.astype(result_dtype) else: result_dtype = result.dtype # Crear perfil de salida profile = dataset.profile.copy() # Determinar valor nodata para el perfil if ignore_nodata: profile_nodata = dataset.nodata else: profile_nodata = dataset.nodata profile.update({ 'dtype': result_dtype, 'nodata': profile_nodata }) # Guardar en archivo si se especifica output_path if output_path: return self._write_raster_to_file(result, profile, output_path) else: return self._create_temp_raster(result, profile)