XYZ to GeoTIFF converter

xyz_to_geotiff.py

# Convert an XYZ file to a GeoTIFF file
#
# Usage
# python xyz_to_geotiff.py input.xyz output.tiff
#
# With custom resolution and different coordinate system
#python xyz_to_geotiff.py input.xyz output.tiff --resolution 0.001 --epsg 3857

import numpy as np
import pandas as pd
from osgeo import gdal, osr
import os

NODATA = 255

def xyz_to_geotiff(xyz_file, output_tiff, resolution=None, epsg=4326):
    """
    Convert XYZ file to GeoTIFF

    Parameters:
    -----------
    xyz_file : str
        Path to input XYZ file
    output_tiff : str
        Path to output GeoTIFF file
    resolution : float, optional
        Resolution of output raster. If None, it's calculated from data
    epsg : int, optional
        EPSG code for the coordinate system (default: 4326, WGS84)
    """
    print(f"Reading XYZ file: {xyz_file}")

    # Read XYZ file (assuming space or tab separated values)
    try:
        # First try comma-separated
        df = pd.read_csv(xyz_file)
        if len(df.columns) < 3:
            # If that fails, try whitespace-separated
            df = pd.read_csv(xyz_file, delim_whitespace=True, header=None)
            if len(df.columns) >= 3:
                df.columns = ['x', 'y', 'z'] + [f'col_{i}' for i in range(3, len(df.columns))]
            else:
                raise ValueError("Could not parse XYZ file - expected at least 3 columns")
    except Exception as e:
        # Try whitespace-separated with no header
        try:
            df = pd.read_csv(xyz_file, delim_whitespace=True, header=None)
            if len(df.columns) >= 3:
                df.columns = ['x', 'y', 'z'] + [f'col_{i}' for i in range(3, len(df.columns))]
            else:
                raise ValueError("Could not parse XYZ file - expected at least 3 columns")
        except Exception as inner_e:
            raise ValueError(f"Failed to read XYZ file: {e}, then: {inner_e}")

    # Get x, y, z columns
    x = df['x'].values
    y = df['y'].values
    z = df['z'].values

    # Determine raster dimensions and geotransform
    x_min, x_max = x.min(), x.max()
    y_min, y_max = y.min(), y.max()

    # Calculate resolution if not provided
    if resolution is None:
        # Estimate resolution from average point spacing
        x_unique = np.unique(x)
        y_unique = np.unique(y)

        if len(x_unique) > 1 and len(y_unique) > 1:
            x_res = np.diff(np.sort(x_unique)).mean()
            y_res = np.diff(np.sort(y_unique)).mean()
            resolution = min(x_res, y_res)
        else:
            # Default resolution if we can't calculate
            resolution = 0.001  # ~ 100m at equator if using degrees

        print(f"Calculated resolution: {resolution}")

    # Calculate grid dimensions
    cols = int(np.ceil((x_max - x_min) / resolution))
    rows = int(np.ceil((y_max - y_min) / resolution))

    # Create empty raster, specifying NODATA as a "fill" value
    raster = np.full((rows, cols), NODATA)

    # Fill raster with XYZ values
    for i in range(len(x)):
        # Calculate raster indices
        col = int((x[i] - x_min) / resolution)
        row = int((y_max - y[i]) / resolution)  # Y is inverted in raster

        # Ensure indices are within bounds
        if 0 <= row < rows and 0 <= col < cols:
            raster[row, col] = z[i]

    # Create geotransform (GDAL uses the upper left corner)
    geotransform = (x_min, resolution, 0, y_max, 0, -resolution)

    # Create GeoTIFF
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(
        output_tiff,
        cols,
        rows,
        1,  # Number of bands
        gdal.GDT_Float32,  # Data type
        options=['COMPRESS=LZW']
    )

    # Set geotransform and projection
    out_ds.SetGeoTransform(geotransform)

    # Set projection
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg)
    out_ds.SetProjection(srs.ExportToWkt())

    # Write data
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(raster)
    out_band.SetNoDataValue(np.nan)

    # Close dataset to write to disk
    out_ds = None

    print(f"Created GeoTIFF: {output_tiff}")

    return output_tiff

if __name__ == "__main__":
    import argparse

    parser = argparse.ArgumentParser(description="Convert XYZ file to GeoTIFF")
    parser.add_argument("xyz_file", help="Input XYZ file path")
    parser.add_argument("output_tiff", help="Output GeoTIFF file path")
    parser.add_argument("--resolution", type=float, default=None,
                        help="Output raster resolution (if not specified, calculated from data)")
    parser.add_argument("--epsg", type=int, default=4326,
                        help="EPSG code for coordinate system (default: 4326, WGS84)")

    args = parser.parse_args()

    xyz_to_geotiff(args.xyz_file, args.output_tiff, args.resolution, args.epsg)