使用GDAL读取Landsat8 GeoTIFF时数组全零及坐标获取求助

First, let’s figure out why your raster is returning all zeros, then move to a straightforward way to extract pixel coordinates.

Fixing All-Zero Pixel Values #

Valid Landsat8 bands rarely have all zeros—let’s dig into the details to diagnose the issue:

1. Check Key Band Metadata
   Start by verifying the band’s core properties to rule out misinterpretation of data:

   from osgeo import gdal
   
   gtif = gdal.Open("FileName.TIF")
   band = gtif.GetRasterBand(1)
   
   # Print critical metadata
   print(f"Data Type: {gdal.GetDataTypeName(band.DataType)}")
   print(f"NoData Value: {band.GetNoDataValue()}")
   min_val, max_val = band.ComputeRasterMinMax()
   print(f"Pixel Min: {min_val}, Pixel Max: {max_val}")
   

   • If min_val and max_val are truly zero, double-check you’re opening the correct file (e.g., not a QA band or corrupted scene).
   • If the NoData value is set to zero, your raster might be mostly empty—but Landsat8 scenes usually have valid data in most areas unless fully cloud-covered.

2. Inspect the Array’s Data Type
   Sometimes incorrect dtype casting can make valid values appear as zero. Confirm the array’s type and sample a small subset:

   arr = band.ReadAsArray()
   print(f"Array dtype: {arr.dtype}")
   print(f"Sample pixels (top-left 5x5):\n{arr[0:5, 0:5]}")
   

   Most Landsat8 bands use uint16 (16-bit unsigned integers), so if your dtype is something like int8, that could indicate a problem (though GDAL usually handles this automatically).

3. Verify File Access & Integrity
   Ensure your file path is correct (try an absolute path if unsure) and the TIF isn’t corrupted. You can also open it in a tool like QGIS to visually confirm if there’s actual data present.

Getting Pixel Coordinates Efficiently #

To map each pixel to its geographic coordinates, use GDAL’s geotransform data. Here are two practical approaches:

1. Loop-Based Method (Small Images) #

Great for testing or small rasters—easy to read and modify:

# Extract geotransform parameters (origin, resolution, skew)
gt = gtif.GetGeoTransform()
x_origin = gt[0]
y_origin = gt[3]
x_res = gt[1]
y_res = gt[5]

# Get raster dimensions
rows = gtif.RasterYSize
cols = gtif.RasterXSize

# Traverse every pixel
for row in range(rows):
    for col in range(cols):
        # Calculate geographic coordinates (skew is usually 0 for Landsat)
        x = x_origin + col * x_res + row * gt[2]
        y = y_origin + row * y_res + col * gt[4]
        # Get the pixel value
        pixel_val = arr[row, col]
        # Use the data as needed (e.g., print, store)
        print(f"Row: {row}, Col: {col}, X: {x:.2f}, Y: {y:.2f}, Value: {pixel_val}")

2. Vectorized Method (Large Images) #

For big Landsat8 scenes, use numpy to avoid slow nested loops—this is way faster:

import numpy as np

# Create grids of column and row indices
col_indices, row_indices = np.meshgrid(np.arange(cols), np.arange(rows))

# Calculate coordinates for all pixels at once
x_coords = x_origin + col_indices * x_res + row_indices * gt[2]
y_coords = y_origin + row_indices * y_res + col_indices * gt[4]

# Stack coordinates with pixel values for easy processing
pixel_data_with_coords = np.stack((x_coords, y_coords, arr), axis=-1)

This method handles large rasters efficiently, which is perfect for Landsat8’s 30m resolution scenes.

内容的提问来源于stack exchange，提问作者Pulkit Agarwal