# -*- coding: utf-8 -*-
"""

Sistema de Predicción de Peligro de Incendios Forestales de México

Índice de Verdor Relativo


Extrae la capa "NDVI_HANTS" (Relleno) del archivo NetCDF de LaNOT
  - Verifica que la orientación vertical esté correcta.
  - Genera el archivo "NDVI.tif", SRC 4326

"""

# Carpetas de trabajo
dir_base = "D:/proyectos/incendios/proceso_diario/DR"
dir_nc = "E:/ftp/lanot/bandas_productos_y_compuestos_goes/NDVI"
shp_limite_nacional = f'{dir_base}/in/limite_nacional/limite_nacional_WGS84.shp'
dir_tif = "H:/2025-10-03_MADTEC/tif_NDVI_HANTS_LaNOT"
dir_tmp = dir_base + "/tmp"

# Importación de librerías
import os
import glob
import shutil
import time
import xarray as xr
import rasterio
from rasterio.transform import from_bounds
import numpy as np

from rasterio.mask import mask
from rasterio.warp import calculate_default_transform, reproject, Resampling
import geopandas as gpd


def netcdf_to_geotiff(netcdf_path, nombre_capa, output_tif):    
    """
    Genera el archivo TIFF con orientación vertical INVERTIDA
    """
    
    try:
        print(" Generando NDVI con orientación VERTICAL INVERTIDA...")
        print("=" * 60)
        
        # 1. Abrir archivo NetCDF
        ds = xr.open_dataset(netcdf_path)
        
        # 2. Extraer datos NDVI
        ndvi_var = ds[nombre_capa]
        ndvi_data = ndvi_var.values
        fill_value = ndvi_var.attrs.get('_FillValue', -9999)
        
        print(f" Datos NDVI - Forma: {ndvi_data.shape}")
        
        # 3. Reducir a 2D
        if len(ndvi_data.shape) == 4:
            ndvi_2d = ndvi_data[0, 0, :, :]
        elif len(ndvi_data.shape) == 3:
            ndvi_2d = ndvi_data[0, :, :]
        else:
            ndvi_2d = ndvi_data
        
        # 4. Obtener coordenadas
        if 'lat' in ds.variables and 'lon' in ds.variables:
            lat = ds['lat'].values
            lon = ds['lon'].values
        elif 'y' in ds.variables and 'x' in ds.variables:
            lat = ds['y'].values
            lon = ds['x'].values
        else:
            # Buscar en dimensiones
            for dim in ndvi_var.dims:
                if dim in ds.variables:
                    if 'lat' in dim or 'y' in dim:
                        lat = ds[dim].values
                    elif 'lon' in dim or 'x' in dim:
                        lon = ds[dim].values
        
        print(f" Coordenadas originales:")
        print(f"   Lat: [{lat[0]:.3f} ... {lat[-1]:.3f}]")
        print(f"   Lon: [{lon[0]:.3f} ... {lon[-1]:.3f}]")
        
        # 5. MANTENER ORIENTACIÓN ORIGINAL (INVERTIDA)
        print("\n MANTENIENDO ORIENTACIÓN INVERTIDA")
        print("   No se aplica np.flipud() - los datos se mantienen como vienen")
        ndvi_invertido = ndvi_2d  # ← NO voltear, mantener invertido
        lat_invertido = lat       # ← NO voltear coordenadas
        
        print(f" Coordenadas invertidas (sin cambios):")
        print(f"   Lat: [{lat_invertido[0]:.3f} ... {lat_invertido[-1]:.3f}]")
        
        # 6. Calcular transformación con coordenadas INVERTIDAS
        lon_min, lon_max = lon.min(), lon.max()
        lat_min, lat_max = lat_invertido.min(), lat_invertido.max()
        
        transform = from_bounds(lon_min, lat_min, lon_max, lat_max, 
                              len(lon), len(lat_invertido))
        
        # 7. Guardar con orientación INVERTIDA
        metadata = {
            'driver': 'GTiff',
            'dtype': ndvi_invertido.dtype,
            'nodata': fill_value,
            'width': len(lon),
            'height': len(lat_invertido),
            'count': 1,
            'crs': 'EPSG:4326',
            'transform': transform,
            'compress': 'lzw'
        }
        
        print(f" Guardando {output_tif} con ORIENTACIÓN INVERTIDA...")
        with rasterio.open(output_tif, 'w', **metadata) as dst:
            dst.write(ndvi_invertido, 1)
            
            # Agregar metadato indicando orientación
            dst.update_tags(
                orientation="vertical_inverted",
                note="Datos mantienen orientación original del NetCDF"
            )
        
        ds.close()
        
        print(" Archivo generado con ORIENTACIÓN VERTICAL INVERTIDA")
        print(" Norte debería estar en la parte INFERIOR del archivo")
        return True
        
    except Exception as e:
        print(f"❌ Error: {e}")
        return False


"""

Recorta la imagen de acuerdo al límite nacional

"""
def recortar_limite(imagen_path, shapefile_path, salida_path):
    # Cargar shapefile y asegurar que esté en 4326
    limite = gpd.read_file(shapefile_path)
    if limite.crs != 'EPSG:4326':
        limite = limite.to_crs('EPSG:4326')
    
    # Abrir imagen original
    with rasterio.open(imagen_path) as src:
        # Si la imagen no está en 4326, reproyectar primero
        if src.crs != 'EPSG:4326':
            # Reproyectar imagen a 4326
            transform, width, height = calculate_default_transform(
                src.crs, 'EPSG:4326', src.width, src.height, *src.bounds
            )
            
            kwargs = src.meta.copy()
            kwargs.update({
                'crs': 'EPSG:4326',
                'transform': transform,
                'width': width,
                'height': height
            })
            
            # Crear imagen temporal reproyectada
            with rasterio.open('temp_4326.tif', 'w', **kwargs) as dst:
                for i in range(1, src.count + 1):
                    reproject(
                        source=rasterio.band(src, i),
                        destination=rasterio.band(dst, i),
                        src_transform=src.transform,
                        src_crs=src.crs,
                        dst_transform=transform,
                        dst_crs='EPSG:4326',
                        resampling=Resampling.bilinear
                    )
            
            # Ahora recortar la imagen reproyectada
            with rasterio.open('temp_4326.tif') as src_4326:
                out_image, out_transform = mask(src_4326, limite.geometry, crop=True, all_touched=True)
                out_meta = src_4326.meta.copy()
        else:
            # La imagen ya está en 4326, recortar directamente
            out_image, out_transform = mask(src, limite.geometry, crop=True, all_touched=True)
            out_meta = src.meta.copy()
    
    # Actualizar metadatos
    out_meta.update({
        "driver": "GTiff",
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    })
    
    # Guardar imagen recortada
    with rasterio.open(salida_path, "w", **out_meta) as dest:
        dest.write(out_image)
    
    print(f"Imagen recortada en EPSG:4326 guardada como: {salida_path}")


#####################################
#
# Proceso Principal
#
#####################################
if __name__ == "__main__":
    fecha = time.strftime("%d/%m/%Y")
    hora = time.strftime("%H:%M:%S")
    print (f"Inicio: {fecha} {hora}\n")

    # Revisa si hay archivos NetCDF sin procesar el último año disponible en el FTP
    anios = glob.glob(f'{dir_nc}/20??')
    dir_ultimo_anio = anios[-1]
    netcdf_files = sorted(glob.glob(f"{dir_ultimo_anio}/*NDVI*_20??????.nc"), reverse=True)
    for netcdf_file in netcdf_files:

        tif_NDVI = dir_tif + "/NDVIR_" + netcdf_file.split("_")[-1].split(".")[0] + ".tif"
        if not os.path.isfile(tif_NDVI):
            print('Generando:', tif_NDVI)
            nombre_capa = "NDVI_HANTS"
            tmp_tif_NDVI = dir_tmp + "/tmp_NDVI.tif"
            target_crs = "EPSG:4326"
                        
            print("#############################")
            print(netcdf_file)
            print(os.path.basename(tmp_tif_NDVI))

            # Extraer NDVI
            netcdf_to_geotiff(netcdf_file, nombre_capa, tmp_tif_NDVI)

            recortar_limite(tmp_tif_NDVI, shp_limite_nacional, tif_NDVI) 

    fecha = time.strftime("%d/%m/%Y")
    hora = time.strftime("%H:%M:%S")
    print (f"\nProceso terminado:{fecha} {hora}")