---
url: /java/gdal.md
---
# GDAL 入门指南

GDAL（Geospatial Data Abstraction Library，地理空间数据抽象库）是处理栅格地理空间数据的开源"瑞士军刀"。它可以读写数百种格式（GeoTIFF、JPEG、PNG、NetCDF、HDF5、ECW、COG 等），支持重投影、波段运算、裁剪、格式转换、金字塔构建等操作，是遥感与 GIS 领域的基石工具。

***

## 1. 安装

### Python 环境

```bash
pip install gdal          # 某些平台可能失败
pip install gdal==$(python -c "from osgeo import gdal; print(gdal.__version__)")
```

更推荐的方式：使用 conda，会自动装好 C++ 运行时依赖。

```bash
conda install -c conda-forge gdal
```

验证安装：

```python
from osgeo import gdal, osr
print(gdal.__version__)   # 例如 3.9.2
```

### 命令行工具

安装 GDAL 后还附带一组常用 CLI 工具：

| 工具 | 用途 |
|------|------|
| `gdalinfo` | 查看数据集元信息（尺寸、CRS、波段数等） |
| `gdal_translate` | 格式转换、简单参数调整 |
| `gdalwarp` | 重投影 / 重采样 |
| `gdalbuildvrt` | 生成虚拟数据集（VRT） |
| `gdaladdo` | 构建概览（overviews/金字塔） |
| `gdaldem` | DEM 高程图派生（坡度、山体阴影等） |
| `gdal_contour` | 从栅格提取等值线 |
| `nearblack` | 去除黑边 |
| `gdal_grid` | 散点转栅格 |

***

## 2. 核心概念

* **Dataset**：一个打开的栅格文件对象，所有操作的入口。
* **Band**：单个波段（如 RGB 中的 R 波段）。通过 `ds.GetRasterBand(i)` 获取。
* **Block**：磁盘上的存储块，`ReadBlock/WriteBlock` 按块 I/O，适合大文件。
* **Overview**：降采样金字塔，加速大范围浏览显示。
* **NoData**：无效值标记，通过 `GetNoDataValue()` 读取。
* **GeoTransform**：仿射变换矩阵，把像素坐标映射到地理坐标：
  ```
  x_geo = gt[0] + x_pix * gt[1] + y_pix * gt[2]
  y_geo = gt[3] + x_pix * gt[4] + y_pix * gt[5]
  ```
* **SRS / CRS**：空间参考系（坐标系），由 `osr.SpatialReference` 表示，常见 WKT / EPSG 编码 / PROJ 字符串互转。

***

## 3. 最小示例：读写 GeoTIFF

### 3.1 读取

```python
from osgeo import gdal

ds = gdal.Open("dem.tif")
if ds is None:
    raise FileNotFoundError("打不开 dem.tif")

print("尺寸:", ds.RasterXSize, "x", ds.RasterYSize)
print("波段数:", ds.RasterCount)
print("数据类型:", gdal.GetDataTypeName(ds.GetRasterBand(1).DataType))
print("GeoTransform:", ds.GetGeoTransform())
print("投影(WKT):", ds.GetProjection())

band = ds.GetRasterBand(1)
data = band.ReadAsArray()            # NumPy 数组 (y, x)
nodata = band.GetNoDataValue()
stats = band.ComputeStatistics(False)  # (min, max, mean, std)
print("统计:", stats)
```

### 3.2 写入

```python
import numpy as np
from osgeo import gdal, osr

arr = np.random.randint(0, 255, size=(256, 256), dtype=np.uint8)

driver = gdal.GetDriverByName("GTiff")
out = driver.Create("out.tif", 256, 256, 1, gdal.GDT_Byte)
out.SetGeoTransform((440720, 30, 0, 3751320, 0, -30))  # UTM 31N 示例
srs = osr.SpatialReference()
srs.ImportFromEPSG(32631)
out.SetProjection(srs.ExportToWkt())
out.GetRasterBand(1).WriteArray(arr)
out.FlushCache()
del out
```

***

## 4. 常用操作

### 4.1 重投影（gdalwarp）

```bash
gdalwarp -t_srs EPSG:4326 -r bilinear in_utm.tif out_wgs84.tif
```

Python：

```python
from osgeo import gdal

src = gdal.Open("in.tif")
dst = gdal.GetDriverByName("GTiff").CreateCopy("out.tif", src, copyOnly=True)

# 用 WarpOptions 做真正的重采样
warp = gdal.AutoCreateWarpedVRT(src, src.GetGeoTransform(), src.GetProjection(),
                                None, "EPSG:4326", gdal.GRPC_BILINEAR)
gdal.GetDriverByName("GTiff").CreateCopy("out_wgs84.tif", warp, copyOnly=True)
```

### 4.2 裁剪子区域

```python
sub = ds.GetRasterBand(1).ReadAsArray(xoff=100, yoff=100, xsize=500, ysize=500)
```

或整幅 ROI：

```bash
gdal_translate -projwin 440720 3751320 441220 3750820 in.tif crop.tif
```

### 4.3 格式转换

```bash
gdal_translate -of JPEG -co "QUALITY=90" in.tif out.jpg
gdal_translate -ot Float32 -scale 0 255 0 1 in.tif normalized.tif
```

### 4.4 拼接 / VRT

```bash
gdalbuildvrt mosaic.vrt scene1.tif scene2.tif scene3.tif
gdalwarp -t_srs EPSG:32631 -r near mosaic.vrt mosaic.tif
```

### 4.5 构建金字塔（Overviews）

```bash
gdaladdo -r average in.tif 2 4 8 16 32
```

Python：

```python
ds = gdal.Open("in.tif", gdal.GA_Update)
for b in range(1, ds.RasterCount + 1):
    ds.BuildOverviews('AVERAGE', [2, 4, 8, 16])
```

### 4.6 多波段合成（RGB → 单文件）

```python
r = gdal.Open("red.tif").GetRasterBand(1).ReadAsArray()
g = gdal.Open("green.tif").GetRasterBand(1).ReadAsArray()
b = gdal.Open("blue.tif").GetRasterBand(1).ReadAsArray()

rgb = np.stack([r, g, b], axis=0)   # (3, h, w)
driver = gdal.GetDriverByName("GTiff")
out = driver.Create("rgb.tif", rgb.shape[1], rgb.shape[2], 3, gdal.GDT_UInt16)
for i, band_arr in enumerate(rgb, start=1):
    out.GetRasterBand(i).WriteArray(band_arr)
out.SetGeoTransform(gdal.Open("red.tif").GetGeoTransform())
out.SetProjection(gdal.Open("red.tif").GetProjection())
out.FlushCache()
```

### 4.7 内存 Dataset（无需落盘）

```python
mem_drv = gdal.GetDriverByName("MEM")
mem_ds = mem_drv.Create("", width, height, bands, gdal.GDT_Float32)
mem_ds.GetRasterBand(1).WriteArray(arr)
# 可直接作为其它函数的输入
```

***

## 5. 坐标系统（OSR）

```python
from osgeo import osr

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)                # WGS84 经纬度
print(srs.ExportToWkt())
print(srs.ExportToProj4())              # +proj=longlat +datum=WGS84 ...

# 任意两种 SRS 互转
srs2 = osr.SpatialReference()
srs2.ImportFromEPSG(32631)              # UTM 31N

transformer = osr.CoordinateTransformation(srs, srs2)
lon, lat = 120.0, 30.0
x, y, _ = transformer.TransformPoint(lon, lat)
print(x, y)                             # UTM 坐标
```

> GDAL 3.x 起内部使用 OGR 的 CRS 对象，`ImportFromEPSG` 已能正确处理多数情况；老代码里的 `SetWellKnownGeogCS("WGS84")` 仍可用但推荐 EPSG 方式。

***

## 6. 性能技巧

1. **分块读取**：大文件避免一次性 `ReadAsArray()`，改用 `ReadBlock(blockX, blockY)` 或 `RasterIO` 带窗口读取。
2. **`RASTERIO_FLAGS`**：`gdal.GRF_ReadOnly`、`gdal.GRF_RastNoData` 等标志位控制行为。
3. **驱动选项（Creation Options）**：写 GeoTIFF 时常用 `-co COMPRESS=LZW -co TILED=YES -co BLOCKSIZE=256`。
4. **缓存**：`gdal.SetCacheUnlimited()` 可放开内存缓存限制（注意机器内存）。
5. **多线程**：`gdalwarp` 自带 `-wo NUM_THREADS=8`；Python 侧可配合 `concurrent.futures` 对分块并行。
6. **Virtual Raster (VRT)**：不落地即可组合多个文件，I/O 按需读取。

***

## 7. 调试与排错

* `gdalinfo -stats file.tif`：快速确认文件能否被识别、是否有 NoData。
* `GDAL_DEBUGON=1`：环境变量打印详细日志。
* `CPL_CURL_VERBOSE=1`：调试网络数据源（HTTP/S3/WMS）。
* 报错 `No such file or directory` 多半是驱动未注册，先 `gdal.AllRegister()`。
* Windows 上 DLL 找不到：确保用 conda 安装的 gdal，或将 `bin` 目录加入 PATH。

***

## 8. 进一步学习

* 官方文档：https://gdal.org/ （API 参考、教程、CLI 手册）
* 官方 Python 教程：https://gdal.org/en/latest/programs/gdal_utils.html
* 配套矢量库 OGR（现为 GDAL 的一部分）：处理 Shapefile / GeoJSON / GPKG
* 高级用法常结合 `rasterio`（Pythonic 封装）或 `xarray`（多维科学数据）

***

*本指南基于 GDAL 3.x API 编写。不同版本间个别函数签名可能有差异，请以你安装的版本对应文档为准。*