""" tiff_to_wgs84.py ---------------- Georeferences an ungeoreferenced TIFF (e.g. a CAD page-layout export) using the bounding box read from a companion DXF file, then reprojects the result to WGS84 (EPSG:4326). Georreferencia um TIFF sem georreferência (p.ex. exportação de layout de página CAD) usando a caixa delimitadora lida de um arquivo DXF complementar e reprojeta o resultado para WGS84 (EPSG:4326). Usage / Uso: python tiff_to_wgs84.py [options] Options: --out Output filename (default: output_wgs84.tiff) --epsg SIRGAS 2000 UTM zone (default: 31982) 31981 = Zone 21S 31982 = Zone 22S 31983 = Zone 23S 31984 = Zone 24S Requirements / Dependências: pip install rasterio pyproj ezdxf numpy """ import argparse import warnings warnings.filterwarnings("ignore") import numpy as np import ezdxf import rasterio from rasterio.crs import CRS from rasterio.transform import from_bounds from rasterio.warp import calculate_default_transform, reproject, Resampling # ── Command-line arguments ──────────────────────────────────────────────────── # Argumentos de linha de comando def parse_args(): parser = argparse.ArgumentParser( description="Georeference a TIFF using a DXF bounding box and reproject to WGS84." ) # Positional: source TIFF / Arquivo TIFF de entrada (obrigatório) parser.add_argument( "tiff_in", help="Path to the source TIFF file (ungeoreferenced CAD export)." ) # Positional: source DXF / Arquivo DXF de entrada (obrigatório) parser.add_argument( "dxf_in", help="Path to the DXF file whose header bounding box will be used for georeferencing." ) # Optional: output filename / Nome do arquivo de saída parser.add_argument( "--out", default="output_wgs84.tiff", help="Output GeoTIFF filename (default: output_wgs84.tiff)." ) # Optional: EPSG code for source CRS / Código EPSG do CRS de origem parser.add_argument( "--epsg", type=int, default=31982, help="EPSG code for the SIRGAS 2000 UTM zone (default: 31982 = Zone 22S)." ) return parser.parse_args() # ── Content area detection ──────────────────────────────────────────────────── # Detecção da área de conteúdo def find_content_bbox(arr: np.ndarray) -> tuple: """ Detect the bounding box of non-white, non-transparent pixels in the image array. This trims the empty page margins typically found in CAD exports. Detecta a caixa delimitadora dos pixels não brancos e não transparentes no array de imagem. Isso remove as margens de página vazias típicas de exportações CAD. Parameters / Parâmetros: arr : numpy array shaped (bands, rows, cols), dtype uint8 Returns / Retorna: (rmin, rmax, cmin, cmax) — pixel row/column indices of content area """ alpha = arr[3] rgb = arr[:3] # A pixel is "content" if it is not fully white AND not fully transparent. # Um pixel é "conteúdo" se não for totalmente branco E não for totalmente transparente. is_content = ~((rgb[0] > 250) & (rgb[1] > 250) & (rgb[2] > 250)) & (alpha > 0) rows_with = np.any(is_content, axis=1) cols_with = np.any(is_content, axis=0) rmin = int(np.where(rows_with)[0][0]) rmax = int(np.where(rows_with)[0][-1]) cmin = int(np.where(cols_with)[0][0]) cmax = int(np.where(cols_with)[0][-1]) return rmin, rmax, cmin, cmax # ── CRS assignment & reprojection ───────────────────────────────────────────── # Atribuição de CRS e reprojeção def reproject_to_wgs84(content: np.ndarray, xmin, ymin, xmax, ymax, src_epsg: int) -> tuple: """ Assign a SIRGAS 2000 UTM CRS to the cropped content array and reproject it to WGS84 (EPSG:4326). Atribui um CRS UTM SIRGAS 2000 ao array de conteúdo recortado e o reprojeta para WGS84 (EPSG:4326). Returns / Retorna: (reprojected_array, dst_transform, dst_crs) """ bands, ch, cw = content.shape src_crs = CRS.from_epsg(src_epsg) # Map the pixel grid exactly to the DXF real-world extent. # Mapeia a grade de pixels exatamente para a extensão real do DXF. src_transform = from_bounds(xmin, ymin, xmax, ymax, cw, ch) dst_crs = CRS.from_epsg(4326) dst_transform, dst_w, dst_h = calculate_default_transform( src_crs, dst_crs, cw, ch, left=xmin, bottom=ymin, right=xmax, top=ymax, ) dst_data = np.zeros((bands, dst_h, dst_w), dtype=content.dtype) for i in range(bands): reproject( source = content[i], destination = dst_data[i], src_transform = src_transform, src_crs = src_crs, dst_transform = dst_transform, dst_crs = dst_crs, resampling = Resampling.lanczos, ) return dst_data, dst_transform, dst_crs # ── Main ────────────────────────────────────────────────────────────────────── def main(): args = parse_args() # ── 1. Read bounding box from DXF header # ── 1. Lê a caixa delimitadora do cabeçalho do DXF print(f"Loading DXF: {args.dxf_in}") doc = ezdxf.readfile(args.dxf_in) xmin = doc.header["$EXTMIN"][0] ymin = doc.header["$EXTMIN"][1] xmax = doc.header["$EXTMAX"][0] ymax = doc.header["$EXTMAX"][1] print(f" Bounding box : ({xmin:.2f}, {ymin:.2f}) -> ({xmax:.2f}, {ymax:.2f})") print(f" Extent : {xmax - xmin:.1f} m x {ymax - ymin:.1f} m") print(f" Source CRS : EPSG:{args.epsg}") # ── 2. Load TIFF and detect content area # ── 2. Carrega o TIFF e detecta a área de conteúdo print(f"Loading TIFF: {args.tiff_in}") with rasterio.open(args.tiff_in) as src: arr = src.read() # shape: (bands, rows, cols) rmin, rmax, cmin, cmax = find_content_bbox(arr) print(f" Full size : {arr.shape[2]} x {arr.shape[1]} px") print(f" Content area : rows {rmin}-{rmax}, cols {cmin}-{cmax}") print(f" Content size : {cmax - cmin + 1} x {rmax - rmin + 1} px") # Crop to content area, discarding empty page margins. # Recorta para a área de conteúdo, descartando margens de página vazias. content = arr[:, rmin:rmax + 1, cmin:cmax + 1] # ── 3. Assign SIRGAS 2000 CRS and reproject to WGS84 # ── 3. Atribui o CRS SIRGAS 2000 e reprojeta para WGS84 print("Reprojecting to WGS84 ...") dst_data, dst_transform, dst_crs = reproject_to_wgs84( content, xmin, ymin, xmax, ymax, args.epsg ) # ── 4. Write output GeoTIFF (LZW compressed) # ── 4. Grava o GeoTIFF de saída (comprimido com LZW) _, dst_h, dst_w = dst_data.shape with rasterio.open(args.out, "w", driver = "GTiff", dtype = "uint8", width = dst_w, height = dst_h, count = dst_data.shape[0], crs = dst_crs, transform = dst_transform, compress = "lzw") as out_file: out_file.write(dst_data) # ── 5. Print summary # ── 5. Exibe resumo print(f"\nSaved: {args.out}") with rasterio.open(args.out) as v: b = v.bounds print(f" CRS : {v.crs}") print(f" Longitude : {b.left:.6f} -> {b.right:.6f}") print(f" Latitude : {b.bottom:.6f} -> {b.top:.6f}") print(f" Size : {v.width} x {v.height} px") if __name__ == "__main__": main()