GDAL 与栅格数据读写实战

Dataset
  ├── RasterXSize, RasterYSize   影像宽高（像元数）
  ├── GeoTransform               6 元组仿射变换（原点、像元尺寸、旋转）
  ├── Projection                 坐标系 WKT 或 EPSG
  ├── Band 1..N                  各波段，各自有数据类型与 NoData
  └── Metadata                   键值对元数据

import rasterio
with rasterio.open("scene.tif") as src:
    print(src.driver)                 # 驱动名，如 GTiff
    print(src.width, src.height)      # 宽高（像元数）
    print(src.crs)                    # 坐标系
    print(src.transform)              # 仿射变换
    print(src.bounds)                 # 地理范围
    print(src.count)                  # 波段数
    print(src.dtypes)                 # 各波段数据类型
    print(src.nodata)                 # NoData 值
    print(src.res)                    # 像元尺寸 (x, y)

import rasterio
import numpy as np

with rasterio.open("scene.tif") as src:
    band = src.read(1)                          # 读取第 1 个波段，返回 2D 数组
    all_bands = src.read()                      # 读取全部波段，返回 3D 数组
    subset = src.read([1, 2, 3])                # 读取指定波段组合
    masked = src.read(1, masked=True)           # 自动应用 NoData 掩膜
    profile = src.profile                        # 复制元数据用于写新文件

print(band.shape, band.dtype)
print("有效像元比例:", np.count_nonzero(~np.isnan(masked.filled(np.nan))) / band.size)

import rasterio
from rasterio.windows import Window

with rasterio.open("big_scene.tif") as src:
    win = Window(col_off=1000, row_off=2000, width=512, height=512)   # 定义窗口
    data = src.read(1, window=win)                                     # 只读该窗口
    # 也可以用地理坐标构造窗口
    from rasterio.windows import from_bounds
    win_geo = from_bounds(116.0, 39.5, 116.5, 40.0, transform=src.transform)
    patch = src.read(1, window=win_geo)

import rasterio
from rasterio.windows import Window

def iter_blocks(path, block=1024):
    with rasterio.open(path) as src:
        for row in range(0, src.height, block):
            for col in range(0, src.width, block):
                h = min(block, src.height - row)
                w = min(block, src.width - col)
                win = Window(col, row, w, h)
                yield win, src.read(1, window=win)

total = 0
for win, block_data in iter_blocks("big_scene.tif", 1024):
    total += block_data.sum()            # 逐块处理，内存恒定
print("总和:", total)

单块内存 = 块宽 × 块高 × 数据类型字节数 × 波段数
例：1024 × 1024 × 4 字节 × 13 波段 ≈ 52 MB
并行 8 块 → 峰值约 420 MB，加上中间变量要留 2 倍余量

import numpy as np
a = np.array([60000], dtype=np.uint16)
print(a.astype(np.int16))          # 溢出：60000 变成负数 -5536
print(a.astype(np.float32))        # 安全：转为浮点
b = a.astype(np.float32) / 10000.0   # 正确做法：先转浮点再运算，最后按需转回

import rasterio
import numpy as np

with rasterio.open("in.tif") as src:
    profile = src.profile.copy()
    profile.update(dtype="float32", count=1, compress="deflate", tiled=True)
    data = src.read([4, 3, 2]).astype("float32")     # 按 RGB 顺序重排波段
    with rasterio.open("out.tif", "w", **profile) as dst:
        dst.write(data)                              # 写入 3 个波段

gdalwarp -t_srs EPSG:32650 -tr 10 10 -r bilinear -of COG in.tif out.tif   # 重投影到 UTM 50N
gdalwarp -t_srs EPSG:4326 -te 116 39.5 117 40.5 in.tif out.tif            # 指定输出范围
gdalwarp -tap -tr 10 10 in.tif out.tif                                    # 对齐到像元网格

import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling

dst_crs = "EPSG:32650"
with rasterio.open("in.tif") as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)   # 计算目标网格
    profile = src.profile.copy()
    profile.update(crs=dst_crs, transform=transform, width=width, height=height)
    with rasterio.open("out.tif", "w", **profile) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform, src_crs=src.crs,
                dst_transform=transform, dst_crs=dst_crs,
                resampling=Resampling.bilinear)

<VRTDataset rasterXSize="10980" rasterYSize="10980">
  <SRS>EPSG:32650</SRS>
  <GeoTransform>300000, 10, 0, 4500000, 0, -10</GeoTransform>
  <VRTRasterBand dataType="UInt16" band="1">
    <NoDataValue>0</NoDataValue>
    <SimpleSource>
      <SourceFilename>tile_01.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="5490" ySize="5490"/>
      <DstRect xOff="0" yOff="0" xSize="5490" ySize="5490"/>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

gdalbuildvrt mosaic.vrt tile_*.tif                       # 用 VRT 拼接多个瓦片
gdalbuildvrt -resolution average mosaic.vrt tile_*.tif   # 统一分辨率后再拼
gdal_translate -of COG mosaic.vrt mosaic.tif             # 落成实体 COG
gdal_merge.py -o merged.tif -n 0 tile_*.tif              # 直接合并为实体文件

gdal_translate -of COG -co COMPRESS=DEFLATE -co BLOCKSIZE=512 -co OVERVIEWS=IGNORE_EXISTING in.tif out.tif
rio cogeo create in.tif out.tif                          # rasterio 生态的等价工具
rio cogeo validate out.tif                               # 验证是否符合 COG 规范
gdalinfo -json out.tif | python -m json.tool             # 检查块布局与概览

import rasterio                    # rasterio 风格：数组导向
with rasterio.open("in.tif") as src:
    arr = src.read(1)

from osgeo import gdal             # 原生 gdal 风格：句柄导向
ds = gdal.Open("in.tif")
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
ds = None                      # 显式释放，触发刷盘

import dask.array as da
import rioxarray

ds = rioxarray.open_rasterio("big_scene.tif", chunks={"x": 2048, "y": 2048})   # 分块懒加载
ndvi = (ds.sel(band=4) - ds.sel(band=3)) / (ds.sel(band=4) + ds.sel(band=3) + 1e-6)
ndvi = ndvi.compute()          # 触发并行计算

export GDAL_CACHEMAX=512                       # 块缓存 512 MB
export GDAL_NUM_THREADS=ALL_CPUS               # 多线程编解码
export GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR  # 加速云上打开

#!/usr/bin/env bash
set -euo pipefail                              # 遇错即停
mkdir -p out
for f in raw/*.tif; do
  base=$(basename "$f" .tif)
  gdalwarp -t_srs EPSG:32650 -tr 10 10 -tap -r bilinear "$f" "out/${base}_utm.tif"
  gdal_translate -of COG -co COMPRESS=DEFLATE "out/${base}_utm.tif" "out/${base}_cog.tif"
  rm "out/${base}_utm.tif"                     # 清理中间文件
  echo "done: $base"
done