栅格数据处理
第七章 栅格数据处理
7.1 栅格数据基础
7.1.1 栅格数据模型
栅格数据由规则网格的像元（像素）组成：
栅格数据结构
├── 头文件信息
│ ├── 行数/列数
│ ├── 像元大小
│ ├── 地理范围
│ ├── 坐标系
│ └── 无数据值
└── 像元值矩阵
├── 波段1
├── 波段2
└── ...
7.1.2 栅格数据类型
| 数据类型 | 取值范围 | 应用场景 |
|---|---|---|
| Byte | 0-255 | 8位影像 |
| Int16 | -32768 to 32767 | 有符号整数 |
| UInt16 | 0-65535 | 无符号整数 |
| Int32 | ±2.1×10⁹ | 32位整数 |
| Float32 | 浮点数 | DEM、连续数据 |
| Float64 | 双精度浮点 | 高精度计算 |
7.1.3 常见栅格格式
| 格式 | 扩展名 | 特点 |
|---|---|---|
| GeoTIFF | .tif | 最通用，支持压缩 |
| ERDAS IMG | .img | Imagine格式 |
| JPEG2000 | .jp2 | 高压缩比 |
| PNG | .png | 无损压缩 |
| ASCII Grid | .asc | 文本格式 |
| NetCDF | .nc | 多维科学数据 |
| HDF | .hdf | 层级数据 |
7.1.4 PyQGIS栅格操作基础
from qgis.core import QgsRasterLayer, QgsProject
# 加载栅格图层
layer = QgsRasterLayer("/path/to/raster.tif", "raster")
if layer.isValid():
QgsProject.instance().addMapLayer(layer)
# 获取栅格信息
print(f"宽度: {layer.width()}")
print(f"高度: {layer.height()}")
print(f"波段数: {layer.bandCount()}")
print(f"范围: {layer.extent()}")
print(f"CRS: {layer.crs().authid()}")
# 获取像元大小
print(f"X分辨率: {layer.rasterUnitsPerPixelX()}")
print(f"Y分辨率: {layer.rasterUnitsPerPixelY()}")
# 获取数据提供者
provider = layer.dataProvider()
print(f"数据类型: {provider.dataType(1)}") # 波段1的数据类型
7.2 栅格图层属性
7.2.1 图层属性对话框
打开方式：
- 双击图层
- 右键 > 属性
- 图层 > 图层属性
主要选项卡：
| 选项卡 | 功能 |
|---|---|
| 信息 | 元数据和统计 |
| 源 | 数据源设置 |
| 符号 | 渲染方式 |
| 透明度 | 透明度设置 |
| 直方图 | 值分布 |
| 渲染 | 渲染选项 |
| 金字塔 | 概览图层 |
| 元数据 | 描述信息 |
| 图例 | 图例设置 |
7.2.2 栅格信息
# 获取详细统计信息
stats = provider.bandStatistics(1) # 波段1
print(f"最小值: {stats.minimumValue}")
print(f"最大值: {stats.maximumValue}")
print(f"平均值: {stats.mean}")
print(f"标准差: {stats.stdDev}")
# 获取波段描述
for band in range(1, layer.bandCount() + 1):
print(f"波段{band}: {provider.generateBandName(band)}")
7.2.3 无数据值设置
# 设置无数据值
provider.setNoDataValue(1, -9999) # 波段1设置nodata为-9999
# 检查某值是否为nodata
if provider.isNoDataValue(1, value):
print("这是无数据值")
7.3 栅格渲染
7.3.1 单波段灰度
from qgis.core import QgsSingleBandGrayRenderer, QgsContrastEnhancement
# 设置单波段灰度渲染
renderer = QgsSingleBandGrayRenderer(provider, 1) # 波段1
# 设置对比度增强
enhancement = QgsContrastEnhancement(provider.dataType(1))
enhancement.setContrastEnhancementAlgorithm(
QgsContrastEnhancement.StretchToMinimumMaximum
)
enhancement.setMinimumValue(0)
enhancement.setMaximumValue(255)
renderer.setContrastEnhancement(enhancement)
layer.setRenderer(renderer)
layer.triggerRepaint()
7.3.2 单波段伪彩色
from qgis.core import (
QgsSingleBandPseudoColorRenderer,
QgsColorRampShader,
QgsRasterShader
)
# 创建颜色渐变
shader = QgsRasterShader()
colorRamp = QgsColorRampShader()
colorRamp.setColorRampType(QgsColorRampShader.Interpolated)
# 定义颜色节点
colorRamp.setColorRampItemList([
QgsColorRampShader.ColorRampItem(0, QColor(0, 0, 255), "低"),
QgsColorRampShader.ColorRampItem(50, QColor(0, 255, 0), "中"),
QgsColorRampShader.ColorRampItem(100, QColor(255, 0, 0), "高")
])
shader.setRasterShaderFunction(colorRamp)
# 应用渲染器
renderer = QgsSingleBandPseudoColorRenderer(provider, 1, shader)
layer.setRenderer(renderer)
layer.triggerRepaint()
7.3.3 多波段彩色合成
from qgis.core import QgsMultiBandColorRenderer
# RGB合成（假设波段4=红，波段3=绿，波段2=蓝）
renderer = QgsMultiBandColorRenderer(provider, 4, 3, 2)
# 设置各波段对比度
for band in [4, 3, 2]:
enhancement = QgsContrastEnhancement(provider.dataType(band))
enhancement.setContrastEnhancementAlgorithm(
QgsContrastEnhancement.StretchToMinimumMaximum
)
if band == 4:
renderer.setRedContrastEnhancement(enhancement)
elif band == 3:
renderer.setGreenContrastEnhancement(enhancement)
else:
renderer.setBlueContrastEnhancement(enhancement)
layer.setRenderer(renderer)
layer.triggerRepaint()
7.3.4 调色板渲染
from qgis.core import QgsPalettedRasterRenderer
# 获取唯一值
classes = QgsPalettedRasterRenderer.classDataFromRaster(provider, 1)
# 创建渲染器
renderer = QgsPalettedRasterRenderer(provider, 1, classes)
layer.setRenderer(renderer)
7.3.5 山体阴影
from qgis.core import QgsHillshadeRenderer
# 设置山体阴影渲染
renderer = QgsHillshadeRenderer(provider, 1, 315, 45)
# 参数：提供者，波段，方位角（度），高度角（度）
renderer.setZFactor(1.0) # 垂直放大系数
layer.setRenderer(renderer)
layer.triggerRepaint()
7.4 栅格计算器
7.4.1 GUI栅格计算器
打开方式：
栅格 > 栅格计算器
界面组件：
- 栅格波段列表
- 运算符按钮
- 表达式输入框
- 输出设置
表达式示例：
# NDVI计算
( "landsat@4" - "landsat@3" ) / ( "landsat@4" + "landsat@3" )
# 条件判断
("dem@1" > 1000) * 1 + ("dem@1" <= 1000) * 0
# 数学运算
sqrt("raster@1" * "raster@1" + "raster@2" * "raster@2")
# 多栅格运算
"raster1@1" + "raster2@1" - "raster3@1"
# 逻辑运算
("raster@1" > 50) AND ("raster@1" < 100)
7.4.2 PyQGIS栅格计算
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
# 准备输入栅格
entries = []
layer = QgsProject.instance().mapLayersByName("dem")[0]
entry = QgsRasterCalculatorEntry()
entry.ref = "dem@1"
entry.raster = layer
entry.bandNumber = 1
entries.append(entry)
# 计算表达式
calc = QgsRasterCalculator(
'"dem@1" * 0.001', # 表达式：转换单位
'/path/to/output.tif', # 输出路径
'GTiff', # 格式
layer.extent(), # 范围
layer.width(), # 宽度
layer.height(), # 高度
entries # 输入
)
result = calc.processCalculation()
if result == 0:
print("计算成功")
7.4.3 使用Processing进行栅格计算
import processing
# GDAL栅格计算
result = processing.run("gdal:rastercalculator", {
'INPUT_A': '/path/to/raster.tif',
'BAND_A': 1,
'FORMULA': 'A * 0.001',
'OUTPUT': '/path/to/output.tif'
})
# NDVI计算示例
result = processing.run("gdal:rastercalculator", {
'INPUT_A': layer, # 红波段
'BAND_A': 4,
'INPUT_B': layer, # 近红外波段
'BAND_B': 5,
'FORMULA': '(B-A)/(B+A)',
'OUTPUT': 'memory:'
})
7.5 栅格分析
7.5.1 地形分析
坡度：
result = processing.run("gdal:slope", {
'INPUT': dem_layer,
'BAND': 1,
'SCALE': 1,
'AS_PERCENT': False,
'OUTPUT': '/path/to/slope.tif'
})
坡向：
result = processing.run("gdal:aspect", {
'INPUT': dem_layer,
'BAND': 1,
'TRIG_ANGLE': False,
'ZERO_FLAT': True,
'OUTPUT': '/path/to/aspect.tif'
})
山体阴影：
result = processing.run("gdal:hillshade", {
'INPUT': dem_layer,
'BAND': 1,
'Z_FACTOR': 1,
'SCALE': 1,
'AZIMUTH': 315,
'ALTITUDE': 45,
'OUTPUT': '/path/to/hillshade.tif'
})
粗糙度：
result = processing.run("gdal:roughness", {
'INPUT': dem_layer,
'BAND': 1,
'OUTPUT': '/path/to/roughness.tif'
})
TRI（地形粗糙指数）：
result = processing.run("gdal:tpitopographicpositionindex", {
'INPUT': dem_layer,
'BAND': 1,
'OUTPUT': '/path/to/tpi.tif'
})
7.5.2 栅格统计
分区统计：
result = processing.run("native:zonalstatisticsfb", {
'INPUT': polygon_layer,
'INPUT_RASTER': raster_layer,
'RASTER_BAND': 1,
'COLUMN_PREFIX': 'stat_',
'STATISTICS': [0, 1, 2, 4, 5, 6], # Count, Sum, Mean, Std, Min, Max
'OUTPUT': 'memory:'
})
栅格图层统计：
result = processing.run("native:rasterlayerstatistics", {
'INPUT': raster_layer,
'BAND': 1,
'OUTPUT_HTML_FILE': '/path/to/stats.html'
})
栅格直方图：
# 获取直方图数据
histogram = provider.histogram(1) # 波段1
for i, count in enumerate(histogram.histogramVector):
value = histogram.minimum + i * histogram.binWidth
print(f"值{value}: {count}个像元")
7.5.3 栅格重分类
使用查找表重分类：
result = processing.run("native:reclassifybytable", {
'INPUT_RASTER': raster_layer,
'RASTER_BAND': 1,
'TABLE': [0, 100, 1, 100, 200, 2, 200, 300, 3], # min1,max1,new1,...
'NO_DATA': -9999,
'RANGE_BOUNDARIES': 0, # 0=min<value<=max
'OUTPUT': '/path/to/reclass.tif'
})
使用表达式重分类：
# 通过栅格计算器实现
expression = '''
CASE
WHEN "raster@1" < 100 THEN 1
WHEN "raster@1" < 200 THEN 2
WHEN "raster@1" < 300 THEN 3
ELSE 4
END
'''
7.5.4 等值线提取
result = processing.run("gdal:contour", {
'INPUT': dem_layer,
'BAND': 1,
'INTERVAL': 10, # 等高距
'FIELD_NAME': 'ELEV',
'CREATE_3D': False,
'IGNORE_NODATA': True,
'OUTPUT': '/path/to/contours.shp'
})
# 等值面（填充等值线）
result = processing.run("gdal:contour_polygon", {
'INPUT': dem_layer,
'BAND': 1,
'INTERVAL': 10,
'OUTPUT': '/path/to/contour_polygons.shp'
})
7.6 栅格转换
7.6.1 格式转换
# GDAL转换
result = processing.run("gdal:translate", {
'INPUT': input_raster,
'TARGET_CRS': 'EPSG:4326',
'DATA_TYPE': 0, # 0=Use input, 1=Byte, 5=Float32等
'OUTPUT': '/path/to/output.tif'
})
7.6.2 投影转换
result = processing.run("gdal:warpreproject", {
'INPUT': input_raster,
'SOURCE_CRS': 'EPSG:4326',
'TARGET_CRS': 'EPSG:32650',
'RESAMPLING': 0, # 0=最近邻, 1=双线性, 2=三次卷积
'TARGET_RESOLUTION': 30,
'OUTPUT': '/path/to/reprojected.tif'
})
7.6.3 重采样
# 修改分辨率
result = processing.run("gdal:warpreproject", {
'INPUT': input_raster,
'TARGET_CRS': '',
'RESAMPLING': 1, # 双线性
'TARGET_RESOLUTION': 10, # 新分辨率
'OUTPUT': '/path/to/resampled.tif'
})
7.6.4 裁剪栅格
按范围裁剪：
result = processing.run("gdal:cliprasterbyextent", {
'INPUT': raster_layer,
'PROJWIN': '100,120,30,45', # xmin,xmax,ymin,ymax
'OUTPUT': '/path/to/clipped.tif'
})
按掩膜裁剪：
result = processing.run("gdal:cliprasterbymasklayer", {
'INPUT': raster_layer,
'MASK': polygon_layer,
'CROP_TO_CUTLINE': True,
'KEEP_RESOLUTION': True,
'OUTPUT': '/path/to/masked.tif'
})
7.6.5 镶嵌
result = processing.run("gdal:merge", {
'INPUT': [raster1, raster2, raster3],
'PCT': False,
'SEPARATE': False,
'OUTPUT': '/path/to/mosaic.tif'
})
# 使用VRT虚拟镶嵌
result = processing.run("gdal:buildvirtualraster", {
'INPUT': [raster1, raster2, raster3],
'RESOLUTION': 0, # 0=平均, 1=最高, 2=最低
'SEPARATE': False,
'OUTPUT': '/path/to/mosaic.vrt'
})
7.6.6 栅格转矢量
多边形化：
result = processing.run("gdal:polygonize", {
'INPUT': raster_layer,
'BAND': 1,
'FIELD': 'DN',
'OUTPUT': '/path/to/polygons.shp'
})
等值线：
result = processing.run("gdal:contour", {
'INPUT': dem_layer,
'BAND': 1,
'INTERVAL': 10,
'OUTPUT': '/path/to/contours.shp'
})
7.6.7 矢量转栅格
result = processing.run("gdal:rasterize", {
'INPUT': vector_layer,
'FIELD': 'value_field',
'UNITS': 1, # 1=像素
'WIDTH': 1000,
'HEIGHT': 1000,
'EXTENT': vector_layer.extent(),
'OUTPUT': '/path/to/rasterized.tif'
})
# 或使用像元大小
result = processing.run("gdal:rasterize", {
'INPUT': vector_layer,
'FIELD': 'value_field',
'UNITS': 0, # 0=地理单位
'WIDTH': 10,
'HEIGHT': 10,
'OUTPUT': '/path/to/rasterized.tif'
})
7.7 栅格金字塔
7.7.1 什么是金字塔
栅格金字塔是预先计算的降采样版本，用于快速显示大型栅格：
原始分辨率 (Level 0)
↓ 2倍降采样
Level 1 (1/2分辨率)
↓ 2倍降采样
Level 2 (1/4分辨率)
↓ ...
7.7.2 创建金字塔
GUI方式：
图层属性 > 金字塔 > 构建金字塔
命令行/PyQGIS：
result = processing.run("gdal:overviews", {
'INPUT': raster_layer,
'LEVELS': '2 4 8 16',
'CLEAN': False,
'RESAMPLING': 0, # 0=最近邻, 1=平均, 2=Gauss
'FORMAT': 0 # 0=内部, 1=外部(.ovr)
})
7.7.3 删除金字塔
result = processing.run("gdal:overviews", {
'INPUT': raster_layer,
'CLEAN': True # 删除金字塔
})
7.8 栅格样式
7.8.1 颜色映射表
from qgis.core import QgsColorRampShader, QgsRasterShader
from qgis.PyQt.QtGui import QColor
# 创建颜色表
colorRamp = QgsColorRampShader()
colorRamp.setColorRampType(QgsColorRampShader.Discrete) # 离散或插值
# 设置颜色项
items = [
QgsColorRampShader.ColorRampItem(1, QColor(255, 0, 0), "类型1"),
QgsColorRampShader.ColorRampItem(2, QColor(0, 255, 0), "类型2"),
QgsColorRampShader.ColorRampItem(3, QColor(0, 0, 255), "类型3"),
]
colorRamp.setColorRampItemList(items)
# 应用到图层
shader = QgsRasterShader()
shader.setRasterShaderFunction(colorRamp)
renderer = QgsSingleBandPseudoColorRenderer(provider, 1, shader)
layer.setRenderer(renderer)
7.8.2 透明度设置
# 设置整体透明度
layer.renderer().setOpacity(0.5)
# 设置无数据值透明
renderer = layer.renderer()
renderer.setNodataColor(QColor(0, 0, 0, 0))
# 设置特定值透明
from qgis.core import QgsRasterTransparency
transparency = QgsRasterTransparency()
# 设置单值透明
pixelList = [
QgsRasterTransparency.TransparentSingleValuePixel(0, 0, 100) # 值,透明度
]
transparency.setTransparentSingleValuePixelList(pixelList)
renderer.setRasterTransparency(transparency)
7.8.3 混合模式
from qgis.PyQt.QtGui import QPainter
# 设置图层混合模式
layer.setBlendMode(QPainter.CompositionMode_Multiply)
# 可用混合模式:
# CompositionMode_SourceOver (正常)
# CompositionMode_Multiply (正片叠底)
# CompositionMode_Screen (滤色)
# CompositionMode_Overlay (叠加)
# ...
7.9 栅格数据导出
7.9.1 导出为文件
# 使用Processing导出
result = processing.run("gdal:translate", {
'INPUT': layer,
'OPTIONS': 'COMPRESS=LZW', # 压缩选项
'DATA_TYPE': 5, # Float32
'OUTPUT': '/path/to/output.tif'
})
7.9.2 导出带样式
# 渲染后导出为图像
from qgis.core import QgsMapSettings, QgsMapRendererSequentialJob
settings = QgsMapSettings()
settings.setLayers([layer])
settings.setExtent(layer.extent())
settings.setOutputSize(QSize(1000, 1000))
job = QgsMapRendererSequentialJob(settings)
job.start()
job.waitForFinished()
image = job.renderedImage()
image.save('/path/to/rendered.png')
7.9.3 Cloud Optimized GeoTIFF (COG)
result = processing.run("gdal:translate", {
'INPUT': layer,
'OPTIONS': 'TILED=YES|COMPRESS=DEFLATE|COPY_SRC_OVERVIEWS=YES',
'OUTPUT': '/path/to/output_cog.tif'
})
# 或使用专门的COG驱动
result = processing.run("gdal:translate", {
'INPUT': layer,
'FORMAT': 'COG',
'OUTPUT': '/path/to/output.tif'
})
7.10 遥感图像处理
7.10.1 NDVI计算
# 计算NDVI (归一化植被指数)
# NDVI = (NIR - Red) / (NIR + Red)
result = processing.run("gdal:rastercalculator", {
'INPUT_A': landsat_layer, # 红波段
'BAND_A': 4,
'INPUT_B': landsat_layer, # 近红外波段
'BAND_B': 5,
'FORMULA': '(B-A)/(B+A)',
'OUTPUT': '/path/to/ndvi.tif'
})
7.10.2 影像分类
非监督分类（K-means）：
result = processing.run("grass7:i.cluster", {
'input': multiband_raster,
'signaturefile': '/path/to/signature',
'classes': 5,
'output': '/path/to/classified.tif'
})
7.10.3 直方图匹配
result = processing.run("saga:histogrammatching", {
'GRID': input_raster,
'REFERENCE': reference_raster,
'METHOD': 0,
'RESULT': '/path/to/matched.tif'
})
7.11 小结
本章详细介绍了QGIS中栅格数据处理的各种操作：
关键要点：
- 理解栅格数据模型和数据类型
- 掌握各种栅格渲染方式
- 熟练使用栅格计算器
- 了解常用地形分析工具
- 能够进行栅格格式转换和处理
- 掌握栅格金字塔和性能优化
栅格数据处理是遥感和空间分析的重要组成部分。
上一章：第06章 矢量数据处理
下一章：第08章 坐标参考系统与投影
浙公网安备 33010602011771号