GDAL—瓦片格式栅格数据创建和修改

瓦片栅格数据

要点：瓦片是物理分区存储方式；实际的写入操作是逻辑分区方式；通过新建瓦片栅格的方式修改瓦片设置

import numpy as np
from osgeo import gdal, osr

# ========== 参数设置 ==========
cols = 25000  # 列数
rows = 20000  # 行数
bands = 1
dtype = gdal.GDT_Float32

# 预期文件大小 = 25000 × 20000 × 4 = 2,000,000,000 字节 ≈ 1.86 GB

output_raster = r""

print(f"正在创建 {cols} × {rows} 的栅格（约 1.9 GB）...")

# ========== 1. 创建空文件 ==========
driver = gdal.GetDriverByName("GTiff")
ds = driver.Create(output_raster, cols, rows, bands, dtype, options=[
    "TILED=YES",  # 瓦片存储（分块读取必需,物理分区存储）
    "BLOCKXSIZE=256",  # 瓦片宽度
    "BLOCKYSIZE=256",  # 瓦片高度
    "COMPRESS=LZW"  # 压缩（可减少实际大小）
])

# ========== 2. 设置坐标和投影 ==========
x_min, x_res = 110.0, 0.01
y_max, y_res = 35.0, -0.01
ds.SetGeoTransform((x_min, x_res, 0, y_max, 0, y_res))

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
ds.SetProjection(srs.ExportToWkt())

# ========== 3. 写入数据（分块写入，避免内存爆炸） ==========
band = ds.GetRasterBand(1)
band.SetNoDataValue(np.nan)

# 分块写入：每次写入 1000 行，峰值内存约 25000×1000×4 ≈ 100MB
block_rows = 1000

for start_row in range(0, rows, block_rows):    # 逻辑分块处理
    end_row = min(start_row + block_rows, rows)
    current_rows = end_row - start_row

    # 生成这一块的数据（示例：随机数）
    # 注意：生成 25000×1000 的随机数组需要约 100MB 内存
    zone_id = start_row // block_rows + 1  # 分区编号从 1 开始
    block_data = np.full((current_rows, cols), float(zone_id), dtype=np.float32)

    # 写入这块数据
    band.WriteArray(block_data, 0, start_row)   # 写入的数据；列方向偏移；行方向偏移

    # 释放内存
    del block_data

    print(f"  已写入行 {start_row} - {end_row} ({(end_row / rows) * 100:.1f}%)")

# ========== 4. 关闭文件 ==========
ds = None

print(f"\n 大文件创建完成: {output_raster}")
print(f"文件大小约 1.9 GB，可作为 Level 3 练习数据")

# ========== 5. 修改（新建文件）瓦片格式 ==========
gdal.Translate(
    "output_tiled.tif",
    "input.tif",
    creationOptions=["TILED=YES", "BLOCKXSIZE=128", "BLOCKYSIZE=128", "COMPRESS=LZW"]
)