import xarray as xr import numpy as np import rasterio from rasterio.transform import from_bounds from rasterio.crs import CRS import os import pandas as pd import datetime # (Keep your create_sample_netcdf function if you are still using it for testing) # It *should* create a proper datetime coordinate, but if you're using your # own problematic NetCDF, this part is crucial. def convert_netcdf_to_cog_sequence( netcdf_path, variable_name, output_base_dir, level_index=0, start_time_index=0, end_time_index=None, compression="DEFLATE", tile_size=256, flip_lat=True ): print(f"\nIniciando conversão de '{netcdf_path}' para sequência de COGs...") print(f"Variável: '{variable_name}', Nível (índice): {level_index}") try: ds = xr.open_dataset(netcdf_path) except FileNotFoundError: print(f"Erro: Arquivo NetCDF não encontrado em '{netcdf_path}'") return except Exception as e: print(f"Erro ao abrir o arquivo NetCDF: {e}") return if variable_name not in ds.data_vars: print(f"Erro: Variável '{variable_name}' não encontrada no dataset.") print(f"Variáveis disponíveis: {list(ds.data_vars.keys())}") ds.close() return data_array = ds[variable_name] if not os.path.exists(output_base_dir): os.makedirs(output_base_dir) print(f"Diretório de saída criado: '{output_base_dir}'") time_dim_exists = 'time' in data_array.dims if time_dim_exists: num_times = len(data_array['time']) if end_time_index is None: end_time_index = num_times - 1 if not (0 <= start_time_index < num_times and 0 <= end_time_index < num_times and start_time_index <= end_time_index): print(f"Erro: Índices de tempo inválidos. Início: {start_time_index}, Fim: {end_time_index}, Total: {num_times}.") ds.close() return time_indices_to_process = range(start_time_index, end_time_index + 1) print(f"Processando passos de tempo de {start_time_index} a {end_time_index}.") # --- IMPORTANT DIAGNOSTICS FOR DEBUGGING --- print(f"DEBUG: data_array['time'] dtype: {data_array['time'].dtype}") if 'units' in data_array['time'].attrs: print(f"DEBUG: data_array['time'] units: {data_array['time'].attrs['units']}") else: print("DEBUG: data_array['time'] has NO 'units' attribute.") # ------------------------------------------ else: time_indices_to_process = [0] print("Dataset não possui dimensão de tempo, convertendo diretamente (assumindo um único passo).") for i, time_idx in enumerate(time_indices_to_process): print(f"\n--- Processando passo de tempo (índice global): {time_idx} ---") current_data_slice = data_array.isel(time=time_idx) if time_dim_exists else data_array if 'level' in current_data_slice.dims: if level_index >= len(current_data_slice['level']): print(f"Erro: Índice de nível {level_index} fora dos limites (0 a {len(current_data_slice['level']) - 1}). Pulando este passo de tempo.") continue current_data_slice = current_data_slice.isel(level=level_index) print(f"Convertendo dados para o nível: {current_data_slice['level'].item()}") else: print("Dataset não possui dimensão de nível.") if 'lat' in current_data_slice.dims and 'lon' in current_data_slice.dims: if list(current_data_slice.dims) != ['lat', 'lon']: current_data_slice = current_data_slice.transpose('lat', 'lon') print("Dimensões reordenadas para (lat, lon).") else: print("Aviso: Dimensões espaciais não são 'lat' e 'lon'. Verifique a ordem das dimensões.") if flip_lat and 'lat' in current_data_slice.dims: current_data_slice = current_data_slice.isel(lat=slice(None, None, -1)) print("Flip vertical aplicado na dimensão 'lat'.") height, width = current_data_slice.shape min_lon = float(current_data_slice.lon.min()) max_lon = float(current_data_slice.lon.max()) min_lat = float(current_data_slice.lat.min()) max_lat = float(current_data_slice.lat.max()) transform = from_bounds(min_lon, min_lat, max_lon, max_lat, width, height) crs = CRS.from_epsg(4326) profile = { "driver": "GTiff", "height": height, "width": width, "count": 1, "dtype": current_data_slice.dtype, "crs": crs, "transform": transform, "compress": compression, "tiled": True, "blockxsize": tile_size, "blockysize": tile_size, "BIGTIFF": "IF_NEEDED", "nodata": current_data_slice.attrs.get("_FillValue", None) } profile.update( { "driver": "GTiff", "tiled": True, "blockxsize": tile_size, "blockysize": tile_size, "compress": compression, "predictor": 2 if compression in ["DEFLATE", "LZW"] else 1, } ) # Constrói o nome do arquivo de saída com o timestamp if time_dim_exists: # --- THE CRITICAL FIX FOR YOUR ERROR --- time_value_from_netcdf = data_array['time'][time_idx].values # Scenario A: xarray loaded it as numpy.datetime64 (this is the ideal case) if isinstance(time_value_from_netcdf, np.datetime64): current_timestamp_dt = pd.to_datetime(time_value_from_netcdf) print(f"DEBUG: Time value is numpy.datetime64: {current_timestamp_dt}") # Scenario B: xarray loaded it as an integer/float (your current problem) elif isinstance(time_value_from_netcdf, (int, float, np.integer, np.floating)): # --- YOU MUST ADAPT THIS PART BASED ON YOUR NETCDF'S 'units' ATTRIBUTE --- # Example 1: if 'units' is 'hours since 2015-02-02 00:00:00' # base_time = datetime.datetime(2015, 2, 2, 0, 0) # current_timestamp_dt = base_time + datetime.timedelta(hours=int(time_value_from_netcdf)) # Example 2: if 'units' is 'days since 1900-01-01 00:00:00' # base_time = datetime.datetime(1900, 1, 1, 0, 0) # current_timestamp_dt = base_time + datetime.timedelta(days=int(time_value_from_netcdf)) # Example 3: If the integer is just an *index* (0, 1, 2...) # and your NetCDF is, for instance, hourly data starting from a fixed date # You'd need to know the actual start_date of your data series # For the given traceback, your file might be 'Eta10_C00_2015020200_TP2M.nc' # and if its internal 'time' coordinate are just integers like 0, 3, 6, 9... # representing hours from 2015-02-02 00:00:00 # Assuming 'hours since 2015-02-02 00:00:00' for demonstration # You MUST replace this with the correct base time and unit from your NetCDF's 'units' attribute! # If your 'units' attribute is different, this line below is the one to change. base_time_epoch = datetime.datetime(2015, 2, 2, 0, 0) time_unit = datetime.timedelta(hours=1) # Or days=1, minutes=1, etc. current_timestamp_dt = base_time_epoch + time_unit * int(time_value_from_netcdf) print(f"DEBUG: Time value is int/float. Converted to: {current_timestamp_dt}") else: # Fallback for unexpected types print(f"ERROR: Unexpected type for time coordinate: {type(time_value_from_netcdf)}. Attempting to use as is.") current_timestamp_dt = time_value_from_netcdf # This will likely fail if not datetime-like timestamp_str = current_timestamp_dt.strftime('%Y%m%d%H') output_cog_filename = f"{variable_name}_{timestamp_str}.tif" else: output_cog_filename = f"{variable_name}_single_slice.tif" full_output_cog_path = os.path.join(output_base_dir, output_cog_filename) print(f"Escrevendo COG em: '{full_output_cog_path}'") try: with rasterio.open(full_output_cog_path, "w", **profile) as dst: dst.write(current_data_slice.values, 1) dst.build_overviews([2, 4, 8, 16, 32], resampling=rasterio.enums.Resampling.average) dst.update_tags(ns='rio_overview', resampling='average') dst.update_tags(**current_data_slice.attrs) if time_dim_exists: dst.update_tags(NETCDF_TIME=current_timestamp_dt.isoformat()) print(f"Conversão para COG concluída com sucesso: '{full_output_cog_path}'") except Exception as e: print(f"Erro ao escrever o arquivo COG '{full_output_cog_path}': {e}") ds.close() print("\nProcessamento de sequência de COGs concluído.") # --- Example Usage (remains the same as previous version) --- if __name__ == "__main__": netcdf_input_file = "Eta10_C00_2015020200_TP2M.nc" output_cogs_directory = "output_cogs_sequence" variable_to_convert = "TP2M" # Create a sample NetCDF file (if you don't have your own) # This sample file *should* create a proper datetime coordinate # create_sample_netcdf(netcdf_input_file) # UNCOMMENT THIS IF YOU NEED A SAMPLE FILE # 2. Convert the NetCDF file to a sequence of COGs convert_netcdf_to_cog_sequence( netcdf_input_file, variable_to_convert, output_cogs_directory, level_index=0, start_time_index=0, end_time_index=None, compression="DEFLATE", tile_size=512 )