如何利用Python绑定高效填充不同分辨率栅格的NoData值？

Hey, let's fix this problem for you! Your goal is to fill the -999 NoData values in your 5m raster B using the 25m raster A, and you want an efficient, automated Python-based solution. Let's break down what went wrong with your previous attempts and then jump to the better approach.

What was wrong with your earlier attempts? #

• The maximum(A,B) approach: This likely failed because GDAL might not recognize -999 as NoData by default. Even after resampling A to 5m, the maximum function would treat -999 as a valid value (instead of NoData) and keep it instead of using A's value.
• Syntax error in your second gdal_calc command: You used a single equals sign (B=-999) which is an assignment in Python, but you need a double equals sign (B==-999) for comparison. That's why you hit the SyntaxError.

Efficient Python GDAL Solution (No Intermediate Resampling) #

This approach avoids creating a full 5m resampled copy of A, which saves time and disk space. We'll only sample values from A where B has NoData, making the process much more efficient.

Here's the complete, ready-to-use code:

import gdal
import numpy as np

# Update these paths to match your file locations
input_5m_raster = "MNT_5m.tif"    # Your 5m raster with -999 NoData
input_25m_raster = "BATI_5m.tif"  # Your 25m source raster
output_filled_raster = "filled_5m.tif"

# Open the 5m raster (B)
ds_b = gdal.Open(input_5m_raster, gdal.GA_ReadOnly)
band_b = ds_b.GetRasterBand(1)
arr_b = band_b.ReadAsArray()
gt_b = ds_b.GetGeoTransform()  # Get geographic transformation details
proj = ds_b.GetProjection()    # Match output projection to B

# Open the 25m raster (A)
ds_a = gdal.Open(input_25m_raster, gdal.GA_ReadOnly)
band_a = ds_a.GetRasterBand(1)
gt_a = ds_a.GetGeoTransform()
arr_a = band_a.ReadAsArray()   # Load A's data into a numpy array

# Find all positions in B where the value is -999 (NoData)
nodata_y, nodata_x = np.where(arr_b == -999)

# Convert 5m pixel coordinates to geographic coordinates (using pixel centers)
x_geo = gt_b[0] + nodata_x * gt_b[1] + gt_b[1] / 2
y_geo = gt_b[3] + nodata_y * gt_b[5] + gt_b[5] / 2

# Convert geographic coordinates to pixel coordinates in the 25m raster
a_x_idx = ((x_geo - gt_a[0]) / gt_a[1]).astype(int)
a_y_idx = ((y_geo - gt_a[3]) / gt_a[5]).astype(int)

# Filter out positions that fall outside the bounds of the 25m raster
valid_mask = (a_x_idx >= 0) & (a_x_idx < ds_a.RasterXSize) & \
             (a_y_idx >= 0) & (a_y_idx < ds_a.RasterYSize)

# Apply the mask to get only valid positions for filling
valid_y = nodata_y[valid_mask]
valid_x = nodata_x[valid_mask]
valid_a_y = a_y_idx[valid_mask]
valid_a_x = a_x_idx[valid_mask]

# Fill B's NoData spots with values from A
arr_b[valid_y, valid_x] = arr_a[valid_a_y, valid_a_x]

# Create the output raster
driver = gdal.GetDriverByName("GTiff")
ds_out = driver.Create(
    output_filled_raster,
    ds_b.RasterXSize,
    ds_b.RasterYSize,
    1,
    band_b.DataType
)
ds_out.SetGeoTransform(gt_b)
ds_out.SetProjection(proj)
out_band = ds_out.GetRasterBand(1)
out_band.WriteArray(arr_b)
out_band.SetNoDataValue(-999)  # Preserve the original NoData value

# Clean up GDAL resources to avoid memory leaks
ds_b = None
ds_a = None
ds_out = None

Fixing Your gdal_calc Command (If You Prefer Command Line) #

If you still want to use the command line approach, fix the comparison operator and explicitly handle NoData:

1. First resample A to 5m (use -r nearest for discrete data, -r bilinear for continuous terrain data):

   gdalwarp -tr 5 5 -r nearest BATI_5m.tif BATI_resampled_5m.tif
   

2. Run gdal_calc with the correct comparison syntax and NoData flag:

   gdal_calc.py -A MNT_5m.tif -B BATI_resampled_5m.tif --outfile=out.tif --calc="A*(B==-999)+B*(B!=-999)" --NoDataValue=-999
   

Why the Python Approach Is Better #

• Faster: We only process the NoData regions instead of the entire raster, and skip writing a large intermediate resampled file.
• Automation-ready: Easy to integrate into larger Python workflows or batch scripts.
• Memory-efficient: For extremely large rasters, you could modify this code to process in chunks, but the above works smoothly for most use cases.

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