如何用Python GDAL将无人机处理后数组生成带投影的TIFF

Got it, let's break down exactly how to add proper projection and georeferencing to your processed drone image array when saving it to a TIFF using Python's GDAL. Since your original drone metadata isn't GDAL-readable, we'll build everything from scratch using the parameters you already have: center GPS coordinates, flight height, pixel dimensions, and the north-aligned array.

Step 1: Choose the Right Projection System #

First, pick a projection that fits your drone data needs—two common options:

• WGS84 (EPSG:4326): Geographic coordinate system using latitude/longitude in degrees. Great for global reference, but less ideal for measuring distances/areas since it’s spherical.
• UTM (Universal Transverse Mercator): Projected coordinate system using meters. Perfect for local drone surveys because it preserves distances and areas accurately. Calculate your UTM zone using the center longitude:
  • Zone number = floor((center_longitude + 180) / 6) + 1
  • Northern Hemisphere: Use EPSG code 326{zone}; Southern Hemisphere: Use 327{zone}.

We’ll use UTM in the example below since it’s more practical for drone work.

Step 2: Calculate the GeoTransform #

GDAL uses a 6-element GeoTransform array to map pixel coordinates to geographic coordinates. The format is:
[top_left_x, pixel_width, rotation_x, top_left_y, rotation_y, pixel_height]

Since your array is north-aligned, rotation values (rotation_x and rotation_y) are 0. Here’s how to compute the rest:

1. Ground Sampling Distance (GSD): This is the real-world distance each pixel represents (e.g., 0.1 meters per pixel). You mentioned you already have this value; if not, calculate it with GSD = (flight_height * sensor_pixel_size) / camera_focal_length.
2. Top-left coordinates: Derive from your center GPS:
   • top_left_x = center_lon_utm - (image_width * GSD) / 2
   • top_left_y = center_lat_utm + (image_height * GSD) / 2
3. Pixel height: GDAL uses a negative value here because image y-axis points downward, while geographic coordinates (like UTM) point upward. So pixel_height = -GSD.

Step 3: Full Python Code Example #

Here’s a complete script tying it all together. We’ll assume your processed array is a NumPy array, north-aligned, and you have all required parameters:

import gdal
import osr
import numpy as np

# ----------------------
# Your input parameters
# ----------------------
processed_array = np.random.rand(500, 500)  # Replace with your actual array (height, width) or (height, width, bands)
center_lon = 116.3972  # Example: Beijing longitude
center_lat = 39.9075   # Example: Beijing latitude
GSD = 0.1  # 0.1 meters per pixel (adjust to your actual value)
output_tiff_path = "georeferenced_drone.tif"

# ----------------------
# 1. Convert center GPS to UTM coordinates
# ----------------------
# Create WGS84 spatial reference
wgs84_srs = osr.SpatialReference()
wgs84_srs.ImportFromEPSG(4326)

# Calculate UTM zone and EPSG code
utm_zone = int((center_lon + 180) / 6) + 1
utm_epsg = 32600 + utm_zone if center_lat > 0 else 32700 + utm_zone

# Create UTM spatial reference
utm_srs = osr.SpatialReference()
utm_srs.ImportFromEPSG(utm_epsg)

# Transform center lat/lon to UTM x/y
transformer = osr.CoordinateTransformation(wgs84_srs, utm_srs)
center_x_utm, center_y_utm, _ = transformer.TransformPoint(center_lat, center_lon)

# ----------------------
# 2. Calculate GeoTransform
# ----------------------
image_height, image_width = processed_array.shape[:2]
num_bands = processed_array.shape[2] if len(processed_array.shape) == 3 else 1

# Compute top-left UTM coordinates
top_left_x = center_x_utm - (image_width * GSD) / 2
top_left_y = center_y_utm + (image_height * GSD) / 2

# Build GeoTransform array
geo_transform = [
    top_left_x,  # Top-left x
    GSD,         # Pixel width (east-west)
    0,           # Rotation (x-axis)
    top_left_y,  # Top-left y
    0,           # Rotation (y-axis)
    -GSD         # Pixel height (north-south, negative for GDAL's downward y-axis)
]

# ----------------------
# 3. Create and write TIFF file
# ----------------------
# Get GDAL TIFF driver
driver = gdal.GetDriverByName('GTiff')

# Match GDAL data type to your array's dtype
if processed_array.dtype == np.uint8:
    gdal_dtype = gdal.GDT_Byte
elif processed_array.dtype == np.float32:
    gdal_dtype = gdal.GDT_Float32
elif processed_array.dtype == np.int16:
    gdal_dtype = gdal.GDT_Int16
else:
    gdal_dtype = gdal.GDT_Float32  # Fallback to float32

# Create the TIFF dataset
dataset = driver.Create(
    output_tiff_path,
    image_width,
    image_height,
    num_bands,
    gdal_dtype
)

# Set GeoTransform and projection
dataset.SetGeoTransform(geo_transform)
dataset.SetProjection(utm_srs.ExportToWkt())

# Write array data to TIFF
if num_bands == 1:
    dataset.GetRasterBand(1).WriteArray(processed_array)
else:
    for band_idx in range(num_bands):
        dataset.GetRasterBand(band_idx + 1).WriteArray(processed_array[:, :, band_idx])

# Flush and close to save changes
dataset.FlushCache()
dataset = None

print(f"Georeferenced TIFF saved to {output_tiff_path}")

Key Notes to Avoid Mistakes #

• Data Type Matching: Ensure the GDAL data type matches your NumPy array’s dtype—using the wrong type can cause data loss or corruption.
• Array Dimensions: GDAL expects (height, width) for single-band or (height, width, bands) for multi-band arrays. If your array is ordered (bands, height, width), transpose it first with processed_array = processed_array.transpose(1, 2, 0).
• UTM vs WGS84: If you prefer WGS84 (degrees), skip the UTM conversion, set the projection to EPSG:4326, and use degree-based GSD (e.g., ~0.0000009 degrees per pixel for 0.1m at the equator). But UTM is always better for local drone mapping.
• GSD Accuracy: Double-check your GSD calculation—small errors here can throw off georeferencing significantly.

内容的提问来源于stack exchange，提问作者Cz Gong