Python空间数据处理 课件 第5-9章：矢量数据处理进阶- 遥感影像处理与应用

第五章：矢量数据处理进阶数据生成属性操作空间查询常用处理空间分析空间数据库学习目标通过本章的学习，我们将掌握以下内容：学习目标能够将测量数据转为空间矢量数据。对矢量数据的属性进行增删改查操作。基于

GDAL

进行空间查询。对矢量数据要素进行简单操作（如裁剪、参数融合、投影变换）。基于

GDAL

进行常见的空间分析操作（缓冲区、网络分析、叠置）。基于

GDAL

操作空间数据库。第五章：矢量数据处理进阶2026年03月

04日2/

33表格数据转为矢量数据下面展示了如何将

CSV

表格数据转为

Shapefile

数据。案例中使用了

Pandas

库读取数据，然后基于

GDAL

矢量数据的基本概念生成

Shapefile。主要步骤使用

Pandas

读取

CSV

文件。使用

ogr.GetDriverByName

创建并初始化

DataSource

对象。创建

WGS84空间参考

(osr.SpatialReference)。创建图层

(Layer)

并逐一定义并添加属性字段

(ogr.FieldDefn)。遍历

DataFrame

行，新建要素

(Feature)

并赋值，同时将其几何坐标(Geometry)

以

WKT

格式创建并附加，最终完成图层保存。第五章：矢量数据处理进阶2026年03月

04日3/

33示例：创建数据集对象和图层读取数据后，建立

DataSource、空间环境，进而创建

Layer

并附加字段（表头）结构。fromosgeoimportogr,

osrimportpandasas

pd34 df=pd.read_csv("cities.csv")5#

1.

创建

DataSource

驱动与文件ds=ogr.GetDriverByName('ESRI

Shapefile').CreateDataSource("cities.shp")8#

2.

从

EPSG

代码建立

WGS84

投影srs=osr.SpatialReference()srs.ImportFromEPSG(4326)12#

3.

创建点状图层layer=ds.CreateLayer('Cities',srs,

ogr.wkbPoint)15#

4.

追加属性字段定义name_fd=ogr.FieldDefn('Name',

ogr.OFTString)name_fd.SetWidth(24)layer.CreateField(name_fd)layer.CreateField(ogr.FieldDefn('Population',

ogr.OFTReal))layer.CreateField(ogr.FieldDefn('Country',

ogr.OFTString))第五章：矢量数据处理进阶2026年03月

04日4/

33示例：创建要素并赋值利用

apply

函数将表格中每行信息转为要素，添加到图层中。#

设置Feature的几何属性Geometrypoint=ogr.CreateGeometryFromWkt(f"POINT({item['longitude']}{item['latitude']})")feature.SetGeometry(point)#创建Featurelayer.CreateFeature(feature)#

新建Feature并且给其属性赋值defadd_feature(item):feature=

ogr.Feature(layer.GetLayerDefn())feature.SetField('Name',

item['name'])feature.SetField('Population',

item['population'])feature.SetField('Country',

item['country'])78910111213#axis=1表示将函数作用于DataFrame的行df.apply(add_feature,axis=1)ds.FlushCache()1718 del

ds第五章：矢量数据处理进阶2026年03月

04日5/

33直接使用

ogr2ogr

工具除了编写

Python

代码，还可以直接使用

GDAL

提供的命令行工具

ogr2ogr

将实测的表格数据转为

Shapefile

矢量数据。ogr2ogr-s_srsEPSG:4326-t_srsEPSG:4326

\-ooX_POSSIBLE_NAMES=LONG

\-ooY_POSSIBLE_NAMES=LAT

\-mapFieldTypeString=Real

\-f"ESRIShapefile"

\StationTemperature.shp

StationTemperature.csv命令参数说明-s_srs

/

-t_srs：分别指定原始和目标空间参考系统（EPSG:4326

即

WGS

84

坐标系）。-ooX_POSSIBLE_NAMES=LONG：指定

CSV

文件中可能作为

𝑥

坐标的字段名为

LONG。-mapFieldType

String=Real：将字段类型从字符串映射为实数类型。第五章：矢量数据处理进阶2026年03月

04日6/

33数据生成属性操作空间查询常用处理空间分析空间数据库属性操作简介（CRUD）对于矢量数据中每个图层中要素的属性数据，读者可以简单将其看作一个二维表格。学过数据库的读者，可能会想到关系数据库的概念。目前空间矢量数据的属性数据大部分就是以关系表的形式进行存储的。而对于关系数据库的操作，常用的就是增删改查（Create、Retrieve、Update、Delete，即CRUD）操作。增

(Create) 删

(Delete)如：添加一个省简称字段

Abbr。 如：删除已添加的多余字段。改

(Update)如：给省名称统一加上“市”或“省”的后缀。查

(Retrieve)如：找出高中数量大于

1

万所的省份。第五章：矢量数据处理进阶2026年03月

04日8/

33属性操作之增：添加简称字段为图层增加表示简称的

Abbr

自定义属性段。核心思路是先建字典映射，再利用CreateField，最后循环读写。1#

必须设置

update=True

允许覆盖写操作2ds=ogr.Open(r"Province.shp",

update=True)3layer=

ds.GetLayer()4names_dict

=

{'北京':

'京',

'新疆':

'新',

'宁夏':'宁'}

#

省略56#

添加一个字符串类型的全新字段7field=ogr.FieldDefn('Abbr',

ogr.OFTString)8field.SetWidth(5)9layer.CreateField(field)1011forfeaturein

layer:1213141516name=feature.GetField('NAME')feature.SetField('Abbr',names_dict.get(name,

''))#

!必不可少：覆盖应用修改完毕的要素让其真正刷入生效layer.SetFeature(feature)第五章：矢量数据处理进阶2026年03月

04日9/

33属性操作之删：删除特定字段我们再尝试把该

Abbr

字段删除掉。删除的方法包括两个步骤：第一步，从属性表中找到该字段的位置索引；第二步，删除该字段。#

1.

寻找字段索引的函数defget_field_index_by_name(layer:ogr.Layer,name:

str):defn:ogr.FeatureDefn=

layer.GetLayerDefn()foriin

range(defn.GetFieldCount()):ifname==

defn.GetFieldDefn(i).GetName():return

iraiseValueError(f'{name}not

found')8#

2.

调用

DeleteField

执行删除ds:ogr.DataSource=ogr.Open(r"Province.shp",

update=True)layer:ogr.Layer=ds.GetLayer()12index=get_field_index_by_name(layer,

'Abbr')layer.DeleteField(index)GDAL

提供的

DeleteField()

方法传入的参数必须是要删除字段的索引编号，所以需要手写一个查询器将其映射出来。第五章：矢量数据处理进阶2026年03月

04日10/

33属性操作之改：更新省市名称后缀这里要更新

NAME

字段：给直辖市名称后添加“市”，自治区后添加“自治区”，特别行政区后添加“特别行政区”，其他的添加“省”。1#

填充属性值2forfeaturein

layer:3name:str=

feature.GetField('NAME')45#

逻辑操作：增加相应的后缀标识6ifnamein('北京','天津','重庆','上海'):7name

+=

'市'8elif

name

in

('内蒙古',

'广西',

'宁夏',

'新疆','西藏'):9name

+=

'自治区'10elif

name

in

('香港',

'澳门'):11name

+=

'特别行政区'12else:13name

+=

'省'1415feature.SetField('NAME',name)16layer.SetFeature(feature)完成思路同样是遍历图层中的每一个

Feature

要素，然后通过

SetField()

方法更新属性值。第五章：矢量数据处理进阶2026年03月

04日11/

33属性操作之查：查询方法矢量数据的属性一般都是以关系表进行保存的，

所以用户可以使用关系数据库查询语言

SQL进行数据查询。GDAL

支持部分SQL

查询功能。属性数据查询是对属性表中存储的数据进行自定义搜索的操作，通常可以通过如下两种方式进行：使用

SQL

查询 遍历要素查询用户可以遍历图层

Layer

中包含的所有

Feature

要素，然后读取要素的属性数据运用

Python

逻辑进行筛选过滤得到自己想要的结果。第五章：矢量数据处理进阶2026年03月

04日12/

33示例：基于

SQL

的查询法使用

SQL

查询是最直白的方式，支持聚合与排序子句。ds:ogr.DataSource=

ogr.Open(fn)layer:ogr.Layer=

ds.GetLayer()3#

选择出中学数量大于1万所的省份query:str=f'SELECTNAME,HighSchoolFROM{layer.GetName()}WHEREHighSchool>

10000'selected:ogr.Layer=ds.ExecuteSQL(query)forfeaturein

selected:print(feature.GetField('NAME'))9#

选择出中学数量最多的省份（使用排序与取首个记录）query:str=f'SELECTNAME,HighSchoolFROM{layer.GetName()}ORDERBYHighSchool

DESC'selected:ogr.Layer=ds.ExecuteSQL(query)print(f"中学数量最多的省份：{selected.GetFeature(0).GetField('NAME')}")对于极值的检索，理论上可使用嵌套的

SELECT

与聚合

MAX

实现，但在

GDAL

中直接使用ORDER

BY

倒序选取第一条

GetFeature(0)

是更合理的做法，避免由于解析器嵌套导致的错误失效。图层名称中尽量不要包含中文。第五章：矢量数据处理进阶2026年03月

04日13/

33示例：基于要素遍历的查询法如果不使用

SQL，我们可以使用原生

Python

中自带的

filter()

与

sorted()

对可迭代对象进行处理。#

使用

filter

函数对要素属性进行过滤selected=list(filter(lambdaf:f.GetField('HighSchool')>10000,

layer))forfeaturein

selected:print(feature.GetField('NAME'))5#

使用

sorted

方法对要素进行自定义排序，这里使用逆序selected_sorted=sorted(layer,key=lambdaf:f.GetField('HighSchool'),

reverse=True)print(selected_sorted[0].GetField('NAME'))print(selected_sorted[0].GetField('HighSchool'))使用第二种遍历的方式更加方便调试一些。如果对

SQL

语言较为熟悉，推荐使用

SQL

这种声明式编程的方式。第五章：矢量数据处理进阶2026年03月

04日14/

33数据生成属性操作空间查询常用处理空间分析空间数据库空间查询概述空间查询是根据地物的空间位置进行查询的一种数据检索方式。OGC

简单要素规范定义了空间几何体之间的基本空间关系：Equals

(相等),

Disjoint

(不相交),Intersects

(相交),Touches

(接触),Crosses

(穿过),Within

(在内部),Contains

(包含),Overlaps

(重叠)使用

SQL

语句：使用支持空间查询的

SQL

语句。但这种方式只对特定数据源有效，某些不支持。使用空间过滤：使用

GDAL

提供的

SetSpatialFilter()

方法。但这种方式主要用于选择给定外接多边形框内的地物，不能实现精确或其他类型的空间查询。读取核心

Geometry

对象

(推荐)：读取每个要素包含的

Geometry

，手动筛选。因

GDAL

内的Geometry

对象几乎实现了所有

OGC

定义关系，灵活性最强。第五章：矢量数据处理进阶2026年03月

04日16/

33示例：基于地理范围的空间查询目标：从省的面状数据中找出湖北省，然后遍历城市的点数据看是否落在湖北省境内。#

获取湖北省作为基底lyr_province:ogr.Layer=

ds_province.GetLayer()ft_hubei

=

next(filter(lambda

f:

'湖北'

in

f.GetField('NAME'),

lyr_province))4lyr_city=ds_city.GetLayer()#

使用

Within()

过滤出落在湖北省境内的所有市selected=

filter(lambdaf:

f.GetGeometryRef().Within(ft_hubei.GetGeometryRef()),lyr_city10 )11forcityin

selected:print(city.GetField('name'))通过获取图层包含的要素集合，使用

Python

内置属性函数与

Geometry

的方法对该集合实行精准筛选。第五章：矢量数据处理进阶2026年03月

04日17/

33示例：基于距离的空间查询目标：找出离武汉市最近的三座城市。1 cities=

ds.GetLayer()2#

用

filter

函数找出武汉市city:

ogr.Feature

=

next(filter(lambda

f:

'武汉'

in

f.GetField('name'),

cities))5#调用

ResetReading()方法特别重要，如果不重置，后面对

Feature的遍历会出错！cities.ResetReading()8#

根据每个市到武汉市的距离进行排序selected=sorted(cities,key=lambdaf:

f.GetGeometryRef().Distance(city.GetGeometryRef()))11#

从下标1开始，因为排序后距离自身的

0

在第

0

位应当排除foriinrange(1,

4):print(selected[i].GetField('name'))特别注意：迭代器的指针经过了拨动以后（如 next），需要进行重置ResetReading()

或重新加载以作进一步检索应用。第五章：矢量数据处理进阶2026年03月

04日18/

33数据生成属性操作空间查询常用处理空间分析空间数据库裁剪

(Clip)在空间分析中常使用代表研究范围的多边形边界去裁剪全局数据集。调用

Layer.Clip()

即可完成干预。注意：裁剪双方必须具有相同投影体系！#

获取原始数据层与裁剪基准层in_lyr=

ogr.Open("River.shp").GetLayer()extent_lyr=

ogr.Open("ShannXi.shp").GetLayer()4#

新建作为输出的空壳图层，保持坐标与几何种类一致out_lyr=

out_ds.CreateLayer('Clipped',in_lyr.GetSpatialRef(),

in_lyr.GetGeomType())8#

执行基于图形范围的裁切剥离in_lyr.Clip(extent_lyr,out_lyr)陕西省内河流提取结果第五章：矢量数据处理进阶2026年03月

04日20/

33投影变换

(Reproject)投影变换是指将地理空间数据从一种坐标参考系统转换到另一种坐标参考系统的过程。在进行投影变换时需要明确变换之前数据的大地参考系和地图投影的定义。方法

1：使用命令及高层封装实现（推荐）使用

ogr2ogr

命令行工具简单准确：1 ogr2ogr-t_srs"+proj=aea+lat_1=25..."China_Projected.shp

China.shp在

Python

中，该命令可通过

gdal.VectorTranslate()

执行同等转换：1 srs_def="""+proj=aea+lat_1=25+lat_2=47+lat_0=30+lon_0=105

..."""2 gdal.VectorTranslate(dst_file,src_file,dstSRS=srs_def,

reproject=True)第五章：矢量数据处理进阶2026年03月

04日21/

33投影变换方法

2：基本

API

手工实现需要手工建立

CreateDataSource

与

dst_layer。根据源文件创建目标的每个属性字段定义。创建转换对象：osr.CoordinateTransformation(src_srs,

dst_srs)。对每个循环源文件的

Geometry

应用

Transform(ctx)，赋予目标新变量并最后写入

ExportToWkt

的

prj

文件。在不需要针对单独特定的个别元素实行干预时，请尽量使用封装

的gdal.VectorTranslate()

函数。第五章：矢量数据处理进阶2026年03月

04日22/

33要素融合

(Dissolve)将具有物理共界联系的多要素合并为一个集合体。调用

Geometry.UnionCascaded()执行兼具拓扑检查的复杂合并。融合县级行政区得到省级外边界#

初始化多面集装箱

wkbMultiPolygongeoms=ogr.Geometry(ogr.wkbMultiPolygon)3forfeatin

in_lyr:feat.geometry().CloseRings()

#

确立闭合形态geoms.AddGeometry(feat.geometry())67#

多要素合并，抹平等价拓扑冗杂相交处union=

geoms.UnionCascaded()10out_feat=

ogr.Feature(defn)out_feat.SetGeometry(union)out_lyr.CreateFeature(out_feat)第五章：矢量数据处理进阶2026年03月

04日23/

33数据生成属性操作空间查询常用处理空间分析空间数据库缓冲区分析

(Buffer)给行政边界添加晕线#

建立外包。此处的距离数位单位基准#

必定等价关联于原图层本身定义的大地网格数值buff=

geometry.Buffer(6500)#

遍历对象为每个面层发起

Buffer

操作forfeaturein

in_lyr:geometry=

feature.GetGeometryRef()4567891011out_feat=

ogr.Feature(def_feat)out_feat.SetGeometry(buff)out_lyr.CreateFeature(out_feat)第五章：矢量数据处理进阶2026年03月

04日25/

33叠置分析

(Overlay)合并两个图层并保留所有输入图层的几何和属性信息。叠置分析是一种将多个地理图层进行叠加和综合处理，以揭示不同图层之间的空间关系和相互作用的方法：联合

(Union) 相交

(Intersection)提取共同区域，结果仅包含两图层重叠的特区。差异

(Difference)从A

图层中去除与B

图层重叠交合所在的部分。对称差

(Sym.

Diff.)也就是提取独立部分，保留两图层不重合的互异孤体模块。第五章：矢量数据处理进阶2026年03月

04日26/

33代码实现：叠置分析依靠

GDAL

中

Geometry

已经承载完成底层数算方法便能立即使用。#

相交部分：对于两组图层可以使用循环判定相交提取输出if

geom2.Intersects(geom1):output.append(geom2.Intersection(geom1))4#

其他核心方法演示（单

Geometry

层级）：#

联合求并union=

poly1.Union(poly2)8#

差异擦除asym_diff=

poly1.Difference(poly2)11#

对称相斥sym_diff=

poly1.SymDifference(poly2)第五章：矢量数据处理进阶2026年03月

04日27/

33网络分析

(Network

Analysis)研究网络节点与边的连通性并优化路由。GDAL

提供

osgeo.gnm

支持点群汇聚拓扑线网并求算出可行解。1 fromosgeoimport

gnm2driver=gdal.GetDriverByName('GNMFile')ds=driver.Create('.',0,0,0,gdal.GDT_Unknown,options=["..."])ds.CopyLayer(lyr_pipes,

'pipes')ds.CopyLayer(lyr_wells,

'wells')7#

转换系统并以双向可逆边权执行相近点连接dn=

gnm.CastToGenericNetwork(ds)ret=dn.ConnectPointsByLines(['pipes',

'wells'],1e-5,1,1,

gnm.GNM_EDGE_DIR_BOTH)12#

拓扑成立后即可运用

Dijkstra

求出第40至60号阀井的最短路长result=dn.GetPath(40,60,

gnm.GATDijkstraShortestPath)Dijkstra

求解最短路径第五章：矢量数据处理进阶2026年03月

04日28/

33数据生成属性操作空间查询常用处理空间分析空间数据库空间数据库概念空间数据库（Spatial

Database）不仅仅能处理常规的二维表格数据，更能高效存储检索位置信息并管理带有索引的点、线、面实体集。开源对象关系数据库

PostgreSQL

提供了庞健的

PostGIS

扩展模块。它遵循OGC

规范，提供丰富的空间操作函数，使其成为最为广泛和功能强大的开源空间数据库系统。安装流程（Ubuntu

下）：sudoaptinstallpostgresql

postgissudoservicepostgresql

startsudo-upostgres

psqlalteruserpostgreswithpassword

'YourPassword'第五章：矢量数据处理进阶2026年03月

04日30/

33命令行批量导入：shp2pgsql通过

shp2pgsql

命令可以轻松将本地的

Shapefile

直接传送转换写入数据库并创建索引（如示例引入全世界地理多边形表述）。#注意中间管道串流技术

|连结

PostgreSQL指令接收段落shp2pgsql-s4326-I"World_Continents"our_world.world_continent

\|psql-hlocalhost-dpostgis_in_action-Upostgres

-W命令参数配置解释-s：指定输入的空间投影参考系编号。-I：指定在新建关系表的对应实体处预先建立空间索引，极大增加并发查检索速度表现。第五章：矢量数据处理进阶2026年03月

04日31/

33通过

GDAL

连接数据库读取图层既然我们能够通过

Python

操作单体本地数据文件，我们也一定可借助

GDAL

标准化组件进行跨域异端挂载读取空间数据库实例。它的返回结果使用依然使用Dataset

去操作：1 #

连接固定的格式配置命令环：2 db_server=""3 str_conn=f"PG:host={db_server}dbname={db_name}user={db_user}

password={db_passwd}"4#

利用

Open

发起协议通信，连入远端！conneciton=

ogr.Open(str_conn)7#

如同抽取独立文件夹一般直接加载其中子表name=

"our_world.world_continent"layer=

conneciton.GetLayer(name)11forfeatin

layer:print(feat.GetField('continent'))1415 del

conneciton第五章：矢量数据处理进阶2026年03月

04日32/

33本章小结本章作为重要的空间要素掌控环节，主要知识点总结归纳如下：掌握使用

Python

进行

CSV

等格式化文档转换为矢量对象。对于扩展的

SQL

查询及原生代码循环获取修改要素进行了实战使用（增删改查）。对于

OGR

对象底层空间

OGC

关系统系能够利用

Within,

Distance

等探底测算与利用。熟练使用 GDAL

提供的跨封装对象实现处理的要素模型融合抽离（

Clip、UnionCascaded）。对空间计算进行系统涉猎：建立

Buffer

发散、重叠异差

Overlay

模型抽取，以及实现建立

gnm

拓扑管道网构建出最短连接点通路模型。基本掌握基于

PostGIS

配置导入并接入现代工业生产级别空间数据库操作架构。第五章：矢量数据处理进阶2026年03月

04日33/

33第六章：栅格数据处理进阶遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习学习目标核心学习目标遥感影像多波段叠加和波段运算遥感图像拼接和裁剪栅格数据投影变换和格式转化栅格数据和

NumPy数组之间的转化基于

NumPy的栅格波段运算基于

GDAL

命令行工具和相应

API

进行栅格数据处理第六章：栅格数据处理进阶2026年03月

05日2/

47遥感影像分级处理大多经过了辐射定标和几何校正，具有基本的空间参考信息。原始

DN(Digital

Number)

值需要用户自行进行大气校正。当卫星数据中心接收到观测数据以后，会对数据进行不同程度的处理：Level-1级别 Level-2

级别大部分经过了大气校正。例如常用的地表反射率数据（Surface

Reflectance）就属于二级以上产品。不论下载哪个级别的数据，波段叠加、影像裁剪

和

影像拼接

等常规预处理通常都是必不可少的。结合

GDAL

命令行工具与

Python

可以实现高效批处理。第六章：栅格数据处理进阶2026年03月

05日3/

47Landsat

8

影像处理背景以下案例使用从

EarthExplorer

下载的

Landsat

8

影像，以西安市土地利用监测为例：Landsat

8

影像存档通常按照

WRS-2

(Worldwide

Reference

System)

坐标系分瓦片

(Tile)

存储。每个瓦片使用

Path

和

Row

的组合编号。叠加西安市矢量边界与

WRS-2

网格可以发现，西安市跨越了三张瓦片区域：P127R036，P127R037

以及

P126R036。因此，处理整幅西安市影像通常分为几个关键步骤：波段叠加

->

拼接

->

裁剪。第六章：栅格数据处理进阶2026年03月

05日4/

47西安市在由于

WRS2

坐标系下的位置图示：西安市边界与涉及的

P127R036、P127R037、P126R036。瓦片之间有重叠区域。第六章：栅格数据处理进阶2026年03月

05日5/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习为什么要进行波段叠加

(Band

Stacking)？Landsat

8

Level-1

数据在分发时，包含各个独立波段的文件。波段叠加的目的为了图像可视化或后续处理的便利，我们往往需要把单个波段的图像进行叠加（Band

Stacking），使之合并为一个包含多波段（如蓝、绿、红、近红外等）的单一文件。使用工具GDAL

库自带了大量高级处理脚本，其中

gdal_merge.py

就可以用来进行波段叠加和影像拼接操作。第六章：栅格数据处理进阶2026年03月

05日7/

47执行环境准备与路径获取#

激活对应的环境condaactivate

osgeo3#

获取

gdal_merge.py

的完整路径(Get-Command

gdal_merge.py).Path通常

gdal_merge.py

并不是系统直接可运行的

.exe，而是一个

Python

脚本：寻找脚本路径

(PowerShell) Python

执行语法找到全路径后，

应当通过指定python

来运行它：1 python

[gdal_merge.py全路径]

[参数]注意：在

Windows

中如果直接输入

gdal_merge.py

执行，可能会导致使用PyCharm

等编辑器打开，而不是直接运行脚本。而在

Mac/Linux

则可以直接运行。第六章：栅格数据处理进阶2026年03月

05日8/

47gdal_merge.py

核心参数gdal_merge.py[-oout_filename][-ofout_format][-co

NAME=VALUE]*[-pspixelsize_xpixelsize_y][-tap][-separate][-q][-v]

[-pct][-ul_lrulxulylrxlry][-init

"value[value .]"][-nnodata_value][-a_nodata

output_nodata_value][-otdatatype][-createonly]

input_files这里要实现波段叠加功能，主要关注以下参数：-separate：核心参数！指示多个输入文件作为多波段依次存储在一个文件里。-oout_filename：输出文件名称。input_files：需要叠加的单个波段文件列表（按顺序排列）。第六章：栅格数据处理进阶2026年03月

05日9/

47波段叠加操作演示下面演示将

P127R037

瓦片的第

2、3、4、5

波段叠加：1 condaactivate

osgeo2#

执行波段叠加python

C:/Users/TheOne/Applications/miniconda3/envs/osgeo/Scripts/gdal_merge.py\-separate

\-o127037.tif

\5678910LC08LC08LC08LC08._B2.TIF

\._B3.TIF

\._B4.TIF

\._B5.TIF第六章：栅格数据处理进阶2026年03月

05日10/

47波段叠加操作演示第六章：栅格数据处理进阶2026年03月

05日11/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习影像拼接

(Mosaic)

简介基于

GDAL

的脚本不会对重叠区域进行复杂的“匀色”处理，通常是按输入顺序产生层叠效应（后输入的覆盖前面的）。可以通过改变输入文件的先后顺序来优化缝合处的表现。拼接是指对具有部分重叠区域的多景影像进行处理，合并为覆盖整个研究区的全景影像。注意事项 可用工具gdal_merge.py

脚本gdalwarp

命令行二进制程序第六章：栅格数据处理进阶2026年03月

05日13/

47使用

gdal_merge.py

拼接与波段叠加不同，拼接时不要使用

-separate

参数。这样会让多幅影像并在同一个二维平面和对应波段中。pythongdal_merge.py-oXiAn-202108-Mosaic.tif

\-n0-a_nodata0

\345.\LC08.\LC08.\LC08.126036_.127037_.127036_._T1\126036.tif

\._T1\127037.tif

\._T1\127036.tif关键参数解析：-n

0：输入影像中作为无效值（NoData）忽略的值。-a_nodata0：输出影像结果中的

NoData

值。接着列举所有需要被拼接的输入文件。第六章：栅格数据处理进阶2026年03月

05日14/

47使用

gdalwarp

拼接gdalwarp

是一个用

C++

写的重采样和扭曲工具。不仅支持影像拼接，还具有更丰富的功能。.\LC08.\LC08.\LC08.126036_.127037_.127036_._T1\126036.tif

\._T1\127037.tif

\._T1\127036.tif

\gdalwarp.exe

\-srcnodata0

\-dstnodata0

\4567XiAn-202108-Mosaic.tif关键参数：-srcnodata

0

和

-dstnodata0：分别设置输入与输出的

NoData

值。输入文件列表接着写，最后一个参数是输出文件路径。注意这个语法与gdal_merge.py

有轻微不同，并不使用

-o。第六章：栅格数据处理进阶2026年03月

05日15/

47GDAL

命令行二进制

VS.

Python

脚本在使用

GDAL

时，理清工具种类避免命令报错：命令行程序（二进制）如

gdalwarp,

gdalinfo,

gdal_translate。在

PATH

中可直接调用。Windows

下加不加

.exe

都可以。Mac/Linux

原生无扩展名。Python

脚本如

gdal_merge.py,gdal_calc.py。Mac/Linux

且加了运行权限并在

PATH

中：可直接执行。Windows：建议带上

python

解释器前缀以及完整路径以确保不出错：python

C:\ .\gdal_merge.py

args。第六章：栅格数据处理进阶2026年03月

05日16/

47拼接效果展示最终将这三个被波段叠加过的影像进行接边操作：注：图中右上方块（P126R036）拍摄于

8

月

27

日，其它拍摄于

8

月

2

日，导致拼接后能看到边缘颜色有一点差异。第六章：栅格数据处理进阶2026年03月

05日17/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习影像裁剪场景已有研究区域的

Shapefile

数据，想沿着边界进行掩膜并裁剪影像矩阵。剪裁（Clipping

或

Cropping）可以极大减小不必要的数据量。常用于只分析特定研究区的场合。通常有两种情形：基于矢量边界 基于经纬度包围框只是想截取研究区域的四角坐标（Bounding

Box，又称

MBR）。这更简单直接。裁剪的主要工具依然推荐极为高效的

gdalwarp。第六章：栅格数据处理进阶2026年03月

05日19/

47场景

1：基于矢量数据裁剪gdalwarp-cutline.\XiAn.shp

\-crop_to_cutline

\-srcnodata0-dstnodata0

\.\XiAn-202108-Mosaic.tif

.\XiAn-202108-AOI.tif核心参数：-cutline：后接矢量边界文件路径。-crop_to_cutline：强制输出结果图像的范围（BBox）缩小到能恰好包围这个矢量的长方形区域并用这个矢量形状进行掩膜裁剪。第六章：栅格数据处理进阶2026年03月

05日20/

47场景

1：基于矢量数据裁剪第六章：栅格数据处理进阶2026年03月

05日21/

47场景

2：基于经纬度范围裁剪1 gdalwarp-te_srsEPSG:4326

\2 -te107.6133.67109.8034.78

\-srcnodata0-dstnodata0

\.\XiAn-202108-Mosaic.tif

.\XiAn-202108-MBR.tif核心参数：-te

xmin

ymin

xmaxymax：确定输出范围框。-te_srs

EPSG:4326：指明

-te

中的坐标参考系是

WGS

84（经纬度）。如果本身就是同坐标系则不需要。第六章：栅格数据处理进阶2026年03月

05日22/

47场景

2：基于经纬度范围裁剪第六章：栅格数据处理进阶2026年03月

05日23/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习栅格地图代数像计算各种植被指数之类的操作是非常普遍的栅格计算。−以归一化植被指数

(NDVI)

为例：NDVI=NIR

+Redgdal_calc.py

计算工具它是基于

NumPy

进行底层支持的一个轻便表达式计算脚本，可以进行任意自定义的像素级算术运算：gdal_calc.py -calc="表达式" -outfile="输出"

\[-A

filename]

[--A_band=n] .第六章：栅格数据处理进阶2026年03月

05日25/

47gdal_calc.py

计算

NDVI

示范假设之前提取了XiAn-202108-AOI.ti（f

包含4

个波段），并得知第四波段是近红外，第三波段是红波段：1 pythongdal_calc.py

\2 -A.\XiAn-202108-AOI.tif -A_band=4

\3 -B.\XiAn-202108-AOI.tif -B_band=3

\4 -outfile=XiAn-NDVI.tif

\5 -calc="(A-B)/(A+B+1e-8)"

\6 -NoDataValue=0 -hideNoData -type=Float32-quiet注意点：-A

和

-B

可以是同一个文件。-calc

里的字母必须对应大写的变量标识。分母加入了极其微小的1e-8

有效避免分母出现纯

0

报错。第六章：栅格数据处理进阶2026年03月

05日26/

47NDVI

计算结果展示得到新的

NDVI

单波段浮点数栅格，其中植被旺盛的区域

NDVI

值较高：第六章：栅格数据处理进阶2026年03月

05日27/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习重投影的原因WGS84(EPSG:4326)

——全球通用的经纬度Web Mercator (EPSG:3857)——网页底图常用CGCS 2000 (EPSG:4490) ——中国常用的地理坐标源影像与矢量分析数据或其它数据集的空间参考往往不一致：常见参考系 投影目的基于平面坐标计算面积、长度或为了保持其他图层重合，必须统一所有的空间参考。第六章：栅格数据处理进阶2026年03月

05日29/

47使用

gdalwarp

进行投影变换gdalwarp

命令也天然支持强大的投影重采样。#

将

XiAn-202108-AOI

重投影至中国国家

2000

大地坐标系gdalwarp-t_srsEPSG:4490

\-dstnodata0

\4 .\XiAn-202108-AOI.tif

\5 .\XiAn-202108-AOI-2000.tif核心参数：-t_srs

EPSG:4490：强制指定目标的

Target

Spatial

Reference

System（坐标参考系）。第六章：栅格数据处理进阶2026年03月

05日30/

47查看影像变换后的元数据对重投影前后的

TIF

分别执行

gdalinfo

就可以看到栅格行列号（分辨率像素大小）和顶点坐标已经发生了巨大变化：1 gdalinfo

XiAn-202108-AOI-2000.tif图示：原影像与

CGCS2000

投影影像的几何差异第六章：栅格数据处理进阶2026年03月

05日31/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习栅格格式转换工具某个特定的行业软件只接受他们内置或者专有格式进行处理；而最常见的通用分发格式一般是

GeoTIFF。gdal_translate

工具它是专门用于数据格式转换的

GDAL

二进制命令。不仅支持极为广泛的几十种目标格式转写，还能够执行诸如灰度转换、像素值标度缩放（Scaling）等内部调整：1 gdal_translate-ofENVI

\2 .\XiAn-202108-AOI.tif

\3 .\XiAn-202108-AOI-ENVI.dat第六章：栅格数据处理进阶2026年03月

05日33/

47格式转换参数与结果1 gdal_translate-ofENVI

\2 .\XiAn-202108-AOI.tif

\3 .\XiAn-202108-AOI-ENVI.dat核心参数解析：-ofENVI：指示输出为

ENVI

软件专属格式。执行结果：当前目录下会生成

XiAn-202108-AOI-ENVI.dat（存储像素二进制值）和XiAn-202108-AOI-ENVI.hdr（相关的

ASCII

编码元数据，描述头文件信息）两个文件体。第六章：栅格数据处理进阶2026年03月

05日34/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习从离散数据到栅格表面栅格数据插值旨在通过已知的离散数据点（如气象站观测温度）推算出未知区域的像元值，生成连续的“表面”。常见

GIS

栅格插值类型最近邻插值

(Nearest

Neighbor)反距离加权插值

(IDW)克里金插值

(Kriging)GDAL

提供了专用的

gdal_grid

命令行。第六章：栅格数据处理进阶2026年03月

05日36/

47使用

gdal_grid

的反距离加权插值gdal_grid-ainvdist:power=2.0:smoothing=1.0

\-ooX_POSSIBLE_NAMES=LONG

\-ooY_POSSIBLE_NAMES=LAT

\-zfieldTEMPERATURE

\StationTemperature.csv

GridedTemperature.tif核心参数：-ainvdist:power=2.0：设置插值算法与幂参数（反距离权重）。-oo：打开输入的CSV

文件时，通过列名推测哪些字段当做经纬度坐标（LONG

/LAT）。-zfield：指定输入文件中应当被估算插值的数值列（如温度TEMPERATURE）。第六章：栅格数据处理进阶2026年03月

05日37/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习多维数组与

Dataset

的联系遥感影像的数据结构在本质上和

NumPy

的

𝑛-维数组十分相似。 维度的对应长和宽各占据一个空间维度（行数、列数）。多波段就是第三个叠加层次（通道深度）。如果是时序卫星数据序列，甚至包含时间维度（四维）。GDAL

库使用底层的

C++

对象

gdal.Dataset

管理数据元信息。我们时常需要利用它的

API

将特定像元获取并转换为

Python

中

np.ndarray

的等效多维数组，进而利用SciPy

或

sklearn

快速分析。第六章：栅格数据处理进阶2026年03月

05日39/

47Python

API：读取

TIF

为

NumPy

数组1importnumpyas

np2fromosgeoimportgdal,gdal_array34fn=r"C:\...\XiAn-202108-AOI.tif"56#

方法一：先打开

Dataset，后读取单个波段或全量波段的

Array7ds=

gdal.Open(fn)8im

=

ds.ReadAsArray() # >numpy.ndarray

形状(4,3938,

6677)910#

也可只读取指定序号为

1

的单独波段11im1=

ds.GetRasterBand(1).ReadAsArray()12print(im1.shape) # >numpy.ndarray

形状(3938,

6677)1314#

方法二：如果只是简单的纯

NumPy

获取，可跳过

Dataset

实例化过程15im2=gdal_array.LoadFile(fn)注意点：

读取到

ndarray

后，空间投影（Transform

和

Projection）元信息并没随数组保留在内存中，它只是纯数字。第六章：栅格数据处理进阶2026年03月

05日40/

47Python

API：NumPy

数组转存至

GeoTiff将处理好的数组放回成含有地理信息的

GeoTiff，步骤较为严格：if

prototype:ds=driver.CreateCopy(fname,prototype)else:#

如无原型参照，需重建其行列数范围以及传入所需的

Transform

六参数与

Projectionds=driver.Create(fname,array.shape[1],

array.shape[0],1,

gdal.GDT_Float32)ds.SetGeoTransform(transform)ds.SetProjection(projection)defsave2img(fname,array,driver='GTiff',

prototype=None,transform=None,projection=None,

nodata=None):#

创建对应的写入器扩展引擎driver=gdal.GetDriverByName(driver)567891011121314151617ds.GetRasterBand(1).WriteArray(array)ifnodataisnotNone:ds.GetRasterBand(1).SetNoDataValue(nodata)ds.FlushCache()第六章：栅格数据处理进阶2026年03月

05日41/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习手动撰写

NDVI

阵列操作除了用

gdal_calc.py

自动化，使用

Python

脚本能进行任何高级定制处理：12345678910111213importnumpyasnpfromnumpyimportmafromosgeoimport

gdalgdal.UseExceptions()ds=gdal.Open('XiAn-202108-AOI.tif')nodata=ds.GetRasterBand(1).GetNoDataValue()#

将红光波段与近红外波段载入为浮点数red=ds.GetRasterBand(3).ReadAsArray().astype(float)nir=

ds.GetRasterBand(4).ReadAsArray().astype(float)#

合成一个排除所有

NoData

的掩膜布尔数组14mask

=

np.logical_or(red =nodata,

nir=

nodata)1516#

利用

NumPy

掩膜广播求植被指数，并收束阈值夹点17ndvi=(nir-red)/(nir+red+

1e-8)18ndvi=np.clip(ndvi,-1,

1)19ndvi[mask]=

-9999第六章：栅格数据处理进阶2026年03月

05日43/

47基于矢量边界的多波段批处理结合前述拼接、裁剪在

Python

环境下同样可以使用极为优雅的

API

gdal.Warp()，从而将命令行动作封箱进循环批处理中：1 inputs=

[2 '126036.tif',3 '127036.tif',4 '127037.tif'5 ]boundary=

'XiAn.shp'output='XiAn-202108-WGS84.tif'8#

第一个参数是输出名，后接输入列表文件序列；#

Python

包装将许多零碎参数打包为命名赋值（kwargs）gdal.Warp(output,inputs,cutlineDSName=boundary,cropToCutline=True,srcNodata=0,dstNodata=0,dstSRS='EPSG:4326')使用

Warp()

不仅避免被终端阻塞，还能捕捉并打印异常方便排查。第六章：栅格数据处理进阶2026年03月

05日44/

47遥感数据处理波段叠加影像拼接影像裁剪波段运算投影变换

(Reprojection)格式转换栅格数据插值栅格数据与多维数组基于

NumPy

的栅格算术分析课后练习本章小结命令行基石gdal_merge.py

和

gdalwarp

涵盖了数据格式大一统的操作。