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 GBoutput_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 = 1000for 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_dataprint(f"  已写入行 {start_row} - {end_row} ({(end_row / rows) * 100:.1f}%)")# ========== 4. 关闭文件 ==========
ds = Noneprint(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"]
)