基于Python的GIS-RS多源数据处理(TIF/SHP/NC/...)【20250630】

栅格数据以规则网格(像素)的数值矩阵表达地理现象，每个单元格代表一个属性值(如高程、温度)。例如卫星影像、数字高程模型、温度分布图。存储格式包括ENVI DAT、GeoTIFF、JPEG、PNG、ASCII Grid等等。

矢量数据是通过几何图形(点、线、面)表示地理实体，基于坐标和拓扑关系存储空间信息。例如城市(点)、河流(线)、国家边界(面)。存储格式包括Shapefile、GeoJSON、KML、GML等。

（1）栅格数据优点：

自然现象表达：直接记录连续或分类数据(如植被覆盖、高程)。
计算高效：适合矩阵运算(如地图代数、坡度计算)。
简单结构：数据存储为规则网格，易于编程处理(如Python的NumPy)。
视觉直观：支持平滑色彩过渡(如遥感影像渲染)。

（2）矢量数据优点：

高精度：几何形状由数学公式定义，可无限缩放而不失真。
存储高效：仅记录关键坐标，数据体积小(尤其适合稀疏数据)。
灵活分析：支持拓扑分析(邻接、连通性)、网络分析(路径规划)、属性查询。
易编辑：可单独修改单个要素(如移动一个点、调整边界)。

（3）矢量化（Raster → Vector）

用途：从栅格中提取边界(如从卫星图提取建筑物轮廓)。
挑战：可能引入锯齿或拓扑错误(需人工修正)。

（4）栅格化（Vector → Raster）

用途：将矢量数据转换为分类栅格(如土地利用图)。
挑战：精度损失(依赖分辨率，如细小多边形可能丢失)。

矢量数据是"几何的艺术"，适合精确、离散的地理实体。栅格数据是"像素的科学"，擅长表达连续、渐变的地理现象。以下主要基于Python阐述栅格数据和矢量数据的多种处理方式。

波段按行交叉格式(BIL)、波段按像元交叉格式(BIP)以及波段顺序格式(BSQ)是三种用来为多波段影像组织影像数据的常见方法。BIL、BIP和BSQ本身并不是影像格式，而是用来将影像的实际像素值存储在文件中的方案。

⛄栅格数据

示例数据：共7各波段（Coastal aerosol、Blue、Green、Red、NIR、SWIR1、SWIR2）

zhengzhou_Landsat_20240517.tif

zhengzhou_Landsat_20240517.dat

👀常用Python库

（1）GDAL——GDAL教程API

GDAL(Geospatial Data Abstraction Library)是一个开源的地理空间数据处理库，广泛应用于GIS软件如ArcGIS、QGIS等。它提供了多种格式的读写支持，包括GeoTIFF、GeoJSON等，并具有众多功能，如数据转换、地理配准。

核心特点：

格式支持广泛：支持200+栅格/矢量格式(GeoTIFF/HDF/NetCDF等)
专业级处理：投影变换、重采样、镶嵌、栅格计算等
命令行工具集：gdal_translate/gdalwarp/gdalinfo等
多语言绑定：Python/C++/Java 接口统一

依赖库：

NumPy(数组支持)
PROJ(坐标转换)
GEOS(几何运算)
底层C++库(libgdal)

GDAL常用函数：GDAL库简介及函数说明

①gdal.Open()：打开数据。

②gdal.GetDriverByName()：获取指定名称的驱动程序(driver)，比如GeoTIFF格式对应的驱动程序"GTiff"。

③gdal.Warp()：实现裁剪、重采样、几何校正、转换格式、投影转换、查看处理进度等操作。

④虚拟内存：允许在内存中创建和处理栅格数据，避免磁盘I/O操作，特别适合临时数据处理或需要高性能的场景。以/vsimem/开头的路径来创建内存中的文件。

…

（2）Rasterio——Rasterio教程API

Rasterio是一个专门用于栅格数据读写操作的Python库。它基于GDAL库进行二次封装，更加符合Python风格，适用于处理和分析各种栅格数据，如GeoTIFF、ENVI和HDF5等格式。Rasterio提供了强大的工具，可以方便地读取、写入和操作栅格数据，提高数据处理效率。

核心特点：

Pythonic API：GDAL的现代Python封装
上下文管理器：自动资源清理
NumPy集成：栅格数据直接转为数组
窗口读写：高效处理大文件
地理上下文保持：自动维护变换矩阵和CRS

依赖库：

GDAL(核心依赖)
NumPy
affine(变换矩阵)
snuggs(表达式求值)
cligj(命令行支持)

Rasterio常用函数：

①rasterio.open()：打开数据。

②rasterio.mask()：根据给定的几何掩膜提取栅格数据，影像裁剪。

③rasterio.merge()：将所有读取的影像合并成一幅影像，影像镶嵌。

④MemoryFile类：允许在内存中创建和处理栅格数据，避免磁盘I/O操作，特别适合临时数据处理或需要高性能的场景。

…

（3）rioxarray——rioxarray教程API

rioxarray是一个基于rasterio的xarray扩展库，专门用于处理地理空间数据。它结合了xarray的强大数据结构和rasterio的地理空间数据处理能力，使得用户可以更简单高效地进行地理空间数据的读取、处理和分析。rioxarray提供了以下核心功能：栅格数据读取、地理坐标系处理、数据切片和裁剪、数据重采样、数据掩膜、数据写入等等。

核心特点：

xarray扩展：为栅格数据添加地理标签
维度感知：处理时间序列/多波段数据
惰性计算：Dask集成支持大数据
无缝互操作：与Rasterio/GeoPandas集成
链式操作：方法链式调用风格

依赖库：

xarray(核心)
rasterio(I/O)
pyproj(坐标系统)
dask(并行计算)
cartopy(可视化)

（3）Pillow

Pillow是Python Imaging Library（PIL）的一个分支，它提供了广泛的图像处理功能，包括图像缩放、旋转、剪裁、颜色转换、滤镜效果等，可以方便地对图像进行操作。

核心特点：

基础图像处理：裁剪/旋转/缩放/滤镜
格式转换：JPEG/PNG/TIFF/BMP等互转
轻量级：无复杂依赖
非地理空间：不保留地理信息
直方图操作：图像增强基础

依赖库：

纯 Python实现 (部分C扩展)
无强制外部依赖

（4）h5py

h5py是一个用于与HDF5文件交互的Python库。HDF5（Hierarchical Data Format version 5）是一种用于存储和组织大型数据集的文件格式，h5py结合了HDF5的强大功能和Python的易用性，使得处理大型数据集变得更加方便和高效。

核心特点：

HDF5接口：处理科学数据格式
层级存储：类似文件系统的数据组织
高效压缩：支持多种压缩过滤器
并行I/O：支持MPI并行读写
大数据优化：分块存储/部分读取

依赖库：

HDF5 C库
NumPy
mpi4py(可选，并行支持)

👀`GeoTIFF`数据

（1）GDAL方式

影像读取
影像写入
基于矢量边界的影像裁剪Warp()。参考博客①、参考博客②
虚拟内存读写数据（文件路径前缀添加/vsimem/）
…

from osgeo import gdal# 读取影像
def read_img(filename):# 打开文件dataset = gdal.Open(filename)# 获取栅格的描述信息im_Description = dataset.GetDescription()# 获取影像中的一个波段，波段的编号从1到dataset.RasterCountim_band = dataset.GetRasterBand(1)print(im_band.GetMinimum())print(im_band.GetMaximum())# 栅格矩阵的列数im_width = dataset.RasterXSize# 栅格矩阵的行数 im_height = dataset.RasterYSize# 波段数im_bands = dataset.RasterCount# 仿射矩阵，返回的是六个参数的tupleim_geotrans = dataset.GetGeoTransform()# 地图投影信息 返回的是WKT格式的字符串im_proj = dataset.GetProjection()# 将数据写成数组，对应栅格矩阵；同时也可以设置获取对应像元的邻域(8邻域)的像元值im_data = dataset.ReadAsArray(0, 0, im_width, im_height).astype(float)# 关闭对象，文件datasetdel datasetreturn im_proj, im_geotrans, im_data, im_height, im_width

from osgeo import gdal# 写入影像(TIFF格式或者ENVI的Dat格式)
def write_img(filename, im_proj, im_geotrans, im_data):# 判断栅格数据的数据类型if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# 判读数组维数if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelse:im_bands, (im_height, im_width) = 1, im_data.shape# (1)创建文件GeoTIFF驱动driver = gdal.GetDriverByName("GTiff")options = ['INTERLEAVE=BAND']  # BSQ格式(通常默认)# 'LINE' is an unexpected value for INTERLEAVE creation option of type string-select.value must be PIXEL or BAND# options = ['INTERLEAVE=LINE']  # BIL格式# options = ['INTERLEAVE=PIXEL'] # BIP格式dataset = driver.Create(filename, im_width, im_height, im_bands, datatype, options)'''# (2)创建文件ENVI驱动driver = gdal.GetDriverByName("ENVI")options = ['INTERLEAVE=BSQ']  # Band Sequential(通常默认)# options = ['INTERLEAVE=BIL']  # Band Interleaved by Line# options = ['INTERLEAVE=BIP']  # Band Interleaved by Pixeldataset = driver.Create(filename, im_width, im_height, im_bands, datatype, options)'''dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数dataset.SetProjection(im_proj)  # 写入投影if im_bands == 1:dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据else:for i in range(im_bands):dataset.GetRasterBand(i + 1).WriteArray(im_data[i,:,:])del dataset

import os
import rasterio# 通过设置虚拟内存，将格式转换（NC→TIF）后的数据保存到虚拟内存中，再读取进行裁剪（基于矢量）
# 此处读取数据的过程暂时省略（代码不全），仅仅展示虚拟路径保存和读取的方式。# 设置GDAL的虚拟路径标识符
temp_path = "/vsimem/"
shp_path = r"D:\Study\...\shp\MinimumBound_buffer.shp"
output_path = r"D:\Study\...\Himawari_pre"# 获取数据
SWR = tep_data.variables["SWR"][:]
SWR_arr = np.asarray(SWR)
# 虚拟路径文件
temp_path1 = temp_path + dates[0] + "_SWR_5km.tif"
write_img(temp_path1, proj, geotransform, SWR_arr)
# rasterio.open(temp_path1)访问不到GDAL的虚拟路径文件
dataset1 = gdal.Open(temp_path1)
outpath1 = os.sep.join([out_path, "SWR", dates[0] + "_SWR_5km_clip.tif"])
clip_result1 = gdal.Warp(outpath1,dataset1,format='GTiff',cutlineDSName=shp_path,cropToCutline=True,dstNodata=0
)
clip_result1.FlushCache()  # 先强制写入磁盘
del clip_result1           # 再释放内存
gdal.Unlink(temp_path1)    # 删除虚拟文件

（2）Rasterio方式

影像读取和写入
基于矢量边界的影像裁剪mask()
影像镶嵌merge()
虚拟内存读写数据(MemoryFile类)
…

import rasterio# 读取影像
filename = r"D:\Desktop\data\郑州\pre\ly\zhengzhou_Landsat_20240517_clip.tif"
with rasterio.open(filename) as dataset:print("文件名:", dataset.name)print("宽度:", dataset.width)print("高度:", dataset.height)print("波段数量:", dataset.count)print("四至范围:", dataset.bounds)print("坐标参考信息:", dataset.crs)print("坐标参考信息wkt:", dataset.crs.wkt)print("仿射地理变换参数:", dataset.transform)# 读取第一个波段的数据band1 = dataset.read(1)print("波段1的形状", band1.shape)# 输出结果
文件名: D:/Desktop/data/郑州/pre/ly/zhengzhou_Landsat_20240517_clip.tif
宽度: 4650
高度: 2802
波段数量: 7
四至范围: BoundingBox(left=657088.5475157951, bottom=3792558.9565376947, right=796588.5475157951, top=3876618.9565376947)
坐标参考信息: EPSG:32649
坐标参考信息wkt: PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32649"]]
仿射地理变换参数: | 30.00, 0.00, 657088.55|
| 0.00,-30.00, 3876618.96|
| 0.00, 0.00, 1.00|
波段1的形状 (2802, 4650)

import rasterio
import geopandas as gpd
from rasterio.mask import mask# 基于矢量边界的影像裁剪
def Shp_Clip_Raster(Shp_path, Raster_path, Clip_path):# 打开矢量数据集outline_shp = gpd.read_file(Shp_path)with rasterio.open(Raster_path) as dataset:# 裁剪栅格数据out_image, out_transform = mask(dataset, outline_shp.geometry, crop=True)# 更新元数据信息out_meta = dataset.meta.copy()out_meta.update({"driver": "GTiff","height": out_image.shape[1],"width": out_image.shape[2],"transform": out_transform})with rasterio.open(Clip_path, "w", **out_meta) as dest:dest.write(out_image)if __name__ == '__main__':Shp_path = r"D:\...\郑州市\郑州市prj.shp"Raster_path = r"D:\...\LC08_L2SP_124036_20240517_20240521_02_T1_SR_ly.tif"Clip_path = r"D:\...\clip\zhengzhou_Landsat_20240517_clip.tif"Shp_Clip_Raster(Shp_path, Raster_path, Clip_path)

import glob
import os
import rasterio
from rasterio.merge import merge# 影像镶嵌
def Raster1_Merge_Raster2(Rasters_path, Mosaic_path):# 打开栅格数据mosaic_files = [rasterio.open(path) for path in Rasters_path]mosaic_image, mosaic_transform = merge(mosaic_files)mosaic_meta = mosaic_files[0].meta.copy()mosaic_meta.update({"driver": "GTiff","height": mosaic_image.shape[1],"width": mosaic_image.shape[2],"transform": mosaic_transform})# 保存镶嵌后的栅格数据with rasterio.open(Mosaic_path, "w", **mosaic_meta) as dest:dest.write(mosaic_image)if __name__ == '__main__':Rasters_Folder = r"D:\...\mos\Rasters"      # Rasters文件夹下包含两个镶嵌的TIF文件Mosaic_path = r"D:\...\mos\LC08_L2SP_124036_20240517_layer1_2_Mosaic.tif"Rasters_path = glob.glob(os.sep.join([Rasters_Folder, '*.tif']))Raster1_Merge_Raster2(Rasters_path, Mosaic_path)

import geopandas as gpd
import rasterio
from rasterio.io import MemoryFile
from rasterio.mask import mask# 虚拟内存读写数据
# 在内存中实现掩膜裁剪处理，并计算保存NDVI结果
def ndvi_computer(shp_path, image_path, ndvi_path):# 读取矢量边界clip_shp = gpd.read_file(shp_path)with rasterio.open(image_path) as dataset:# 直接在内存中裁剪with MemoryFile() as memfile:out_image, out_transform = mask(dataset, clip_shp.geometry, crop=True)# 更新元数据信息out_meta = dataset.meta.copy()out_meta.update({"driver": "GTiff","height": out_image.shape[1],"width": out_image.shape[2],"transform": out_transform})with memfile.open(**out_meta) as clipped_ds:clipped_ds.write(out_image)print("裁剪后尺寸:", out_image.shape)# 使用内存栅格进行分析with memfile.open() as ds:ndvi = (ds.read(4) - ds.read(3)) / (ds.read(4) + ds.read(3))# 保存NDVI结果save_ndvi(ndvi, ds.meta, ndvi_path)def save_ndvi(ndvi, image_meta, ndvi_path):# 更新元数据信息ndvi_meta = image_meta.copy()ndvi_meta.update({"driver": "GTiff","count": 1,"dtype": "float32","nodata": 0})# 写入文件with rasterio.open(ndvi_path, "w", **ndvi_meta) as dest:# 写入数据到第1波段dest.write(ndvi, 1)if __name__ == '__main__':shp_path = r"D:\Desktop\data\...\Layer1.shp"image_path = r"D:\Desktop\data\...\Landsat_20240517_clip.tif"ndvi_path = r"D:\Desktop\data\...\Layer1_NDVI.tif"# 裁剪影像并计算NDVI，同时将裁剪后的影像先保存在内存中，再读取ndvi_computer(shp_path, image_path, ndvi_path)

小编基于rasterio库的"mask"和ArcGIS的"Extract by Mask"工具两种方式，对栅格影像进行裁剪：

利用ArcGIS的"Extract by Mask"工具裁剪影像会出现像元偏移的问题，像元值与原始影像存在比较小的差异；利用rasterio库的"mask"裁剪影像，像元值是对应的。
基于rasterio库的"mask"和ArcGIS的"Extract by Mask"工具两种方式对原始影像进行裁剪，分别得到有交集的Layer1和Layer2，统一利用rasterio库的"merge"进行影像镶嵌，同样存在上述问题。

在Rasterio中，dataset.profile和dataset.meta都是用于访问栅格数据集元数据的属性，但它们在使用场景和返回内容上有一些关键区别：rasterio中dataset.profile和dataset.meta的区别：

特性	dataset.profile	dataset.meta
返回类型	字典(可写，可直接用于创建新文件)	字典(只读副本)
主要用途	创建新栅格时的模板	查看元数据
可变性	可修改(修改后会影响数据集)	不可修改
包含内容	核心元数据 + 部分驱动特定选项	仅核心元数据
版本引入	`Rasterio 1.0+`	早期版本即有

rasterio中dataset.meta属性设置及其详细说明：

属性名	数据类型	说明
`driver`	`str`	栅格驱动名称(如`'GTiff'`、`'PNG'`、`'HFA'`等)
`dtype`	`str`	数据类型(如`'uint8'`、`'float32'`等)，对应`NumPy`数据类型
`nodata`	`int`/`float`/`None`	表示无效数据的值(如`0`、`-9999`或`None`表示无`Nodata`值)
`width`	`int`	栅格列数(宽度)
`height`	`int`	栅格行数(高度)
`count`	`int`	波段数量(如`RGB`图像为3)
`crs`	`rasterio.crs.CRS`	坐标参考系统(如`CRS.from_epsg(4326)`表示`WGS84`)
`transform`	`Affine`	仿射变换矩阵

# 示例meta输出
{'driver': 'GTiff','dtype': 'uint16','nodata': None,'width': 1000,'height': 800,'count': 3,'crs': CRS.from_epsg(32650),'transform': Affine(10.0, 0.0, 500000.0,0.0, -10.0, 1200000.0)
}
# 示例profile输出
{'driver': 'GTiff','dtype': 'uint16','nodata': 0.0,'width': 1024,'height': 768,'count': 3,'crs': CRS.from_epsg(4326),'transform': Affine(0.1, 0.0, 30.0,0.0, -0.1, 50.0),'compress': 'lzw',          # 驱动特定选项'tiled': True,              # 驱动特定选项'blockxsize': 256,          # 驱动特定选项'blockysize': 256           # 驱动特定选项
}

使用场景示例：

①创建新文件时：优先使用.profile

②修改参数时：复制.profile后再修改

③仅查看元数据时：使用.meta更安全

④需要驱动特定参数时：必须用.profile

总结来说：需要修改或创建新文件时用.profile，只读访问时用.meta。

👀`HDF`数据

HDF(Hierarchical Data Format)是用于存储和分发科学数据的一种自我描述、多对象文件格式。HDF是由美国国家超级计算应用中心(NCSA)创建的，中文名称为层次型数据格式，以满足不同群体的科学家在不同工程项目领域之需要。HDF主要特征是：

自述性： HDF文件里的每一个数据对象，有关于该数据的综合信息（元数据）。在没有任何外部信息的情况下，应用程序能够解读HDF文件的结构和内容。
通用性：许多数据类型都可以被嵌入在一个HDF文件里。例如，通过使用合适的HDF数据结构，符号、数字和图形数据可以同时存储在一个HDF文件里。
灵活性：HDF允许用户把相关的数据对象组合在一起，放到一个分层结构中，可在数据对象中添加描述和标签。它还允许用户把科学数据放到多个HDF文件里。
扩展性：HDF极易容纳新增加的数据模式，容易与其他标准格式兼容。
跨平台性：HDF是一个与平台无关的文件格式。HDF文件无需任何转换就可以在不同平台上使用。

HDF的最初产生于20世纪80年代，到现在已经具有两个不同的产品。从HDF1到HDF4的各个版本在本质上是一致的，因而HDF4可以兼容早期的版本。HDF5推出于1998年，相较于以前的HDF文件，可以说是一种全新的文件格式。它与HDF4只在概念上一脉相承，而在数据结构的组织上却截然迥异。HDF5的产生与发展反映了HDF在不断适应现代计算机发展和数据处理日益庞大复杂的要求。HDF强大的机制适应了遥感影像的特点，能够有条不紊、完备地保存遥感影像的属性和空间信息数据，同时使查询和提取相关数据也很方便容易。HDF4与HDF5的优缺点对比：

HDF4文件由文件头，数据描述符块和数据元素组成，后两者组成数据对象。数据描述符块由若干描述符组成，它们由标识符、参照符、数据偏移量、数据长度等组成。标识符和参照数组合在一起唯一确定一个数据对象。
HDF4不能存储多于2万个复杂对象，文件大小不能超过2G字节，其数据结构不能完全包含它的内容。随着对象的增多，数据类型也受到限制，库代码过时，过于琐碎，不能有效执行并行I/O，难于运用到线程程序中，HDF5不但能处理更多对象，存储更大的文件，支持并行I/O，线程和具备现代操作系统与应用程序所要求的其他特性。而且数据模型变得更简单，概括性更强。
HDF5格式运用了HDF4和AIO文件的某些关键思想, 比HDF4的自描述性更强, 它由一个超级块(super block)、B树结点(B-tree node)、对象头(object header)、集合(collection)、局部堆(local heaps)和自由空间(free space)组成。

HDF是一种功能强大, 广泛运用于科学领域的文件格式。研究它的组织结构特别是HDF5的组织结构对于处理和管理地理信息系统的海量图形数据和属性数据具有一定的借鉴作用。掌握和运用NCSA提供的API提取影像数据，可以节省时间，提高程序编写效率。因此，HDF将会得到很广泛的应用。而HDF5由于前面所述很多优点，已经可以在未来取代HDF4。

# 小编以批处理夜间灯光数据VNP46A4为例
# VNP46A4.A2012001.h28v04.001.2021125015952.hdf
# VNP46A4.A2012001.h28v05.001.2021125045353.hdfimport os
import h5py
import numpy as np
from osgeo import gdal, osr# 保存影像
def write_img(filename, im_proj, im_geotrans, im_data):# 判断栅格数据的数据类型if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# 判读数组维数if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelse:im_bands, (im_height, im_width) = 1, im_data.shape# 创建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数dataset.SetProjection(im_proj)  # 写入投影if im_bands == 1:dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据else:for i in range(im_bands):dataset.GetRasterBand(i + 1).WriteArray(im_data[i,:,:])del datasetdef hdf_to_tif(input_folder_path, output_path):# 批量获取hdf文件名files_name = []# 获取目标路径下的文件和文件夹的名字列表for path in os.listdir(input_folder_path):if os.path.isfile(os.path.join(input_folder_path, path)):files_name.append(path)for i in range(len(files_name)):input_path = os.sep.join([input_folder_path, files_name[i]])with h5py.File(input_path, "r") as f:'''for key in f.keys():# print(f[key], key, f[key].name, f[key].value)# 因为这里有group对象,它是没有value属性的,故会异常。另外字符串读出来的是字节流，需要解码成字符串# f[key]表示dataset或group对象,f[key].value访问dataset的值,group对象除外。print(f[key], key, f[key].name)print("*"*50)'''# 根据夜间灯光数据hdf格式的数据结构,读取对应的数据层获取数据HDF5_group = f["HDFEOS"]["GRIDS"]["VIIRS_Grid_DNB_2d"]["Data Fields"]# print(HDF5_group.keys())DNB_data = HDF5_group["AllAngle_Composite_Snow_Free"][:]*0.1# 获取hdf文件中对应经纬度变量的信息lon_data = HDF5_group["lon"][:]lat_data = HDF5_group["lat"][:]# 影像的左上角&右下角坐标lonmin, latmax, lonmax, latmin = [lon_data.min(), lat_data.max(), lon_data.max(), lat_data.min()]# 分辨率计算num_lon = len(lon_data)num_lat = len(lat_data)lon_res = (lonmax - lonmin) / (float(num_lon) - 1)lat_res = (latmax - latmin) / (float(num_lat) - 1)# 定义投影proj = osr.SpatialReference()proj.ImportFromEPSG(4326)  # WGS84proj = proj.ExportToWkt()  # 重点，转成wkt格式# 定义六参数，设置影像的显示范围和分辨率# 影像左上角横坐标：geoTransform[0]# 影像左上角纵坐标：geoTransform[3]# 遥感图像的水平空间分辨率为geoTransform[1]# 遥感图像的垂直空间分辨率为geoTransform[5]# 通常geoTransform[5]与geoTransform[1]相等# 如果影像方向没有发生旋转，即上北、下南，则geoTransform[2]与geoTransform[4]为零。geotransform = (lonmin, lon_res, 0.0, latmax, 0.0, -lat_res)# 注意os.path.join和os.sep.join的区别save_name = files_name[i].split(".")save_name = "_".join(save_name[:3])output_file_path = os.path.join(output_path, save_name + ".tif")write_img(output_file_path, proj, geotransform, DNB_data)print("转换成功")if __name__ == '__main__':# hdf文件输入输出路径input_folder_path = r"D:\Desktop\data\hdf\VNP46A4"output_path = r"D:\Desktop\data\hdf\tif"# 批量读取hdf文件，转换为tif文件hdf_to_tif(input_folder_path, output_path)

👀`DAT`数据

ENVI软件中的.dat文件是存储遥感图像原始像元值数据的二进制文件。它本身缺乏解释自身内容的信息，必须依赖一个同名的.hdr头文件来提供解读这些二进制数据所需的元数据(尺寸、波段数、数据类型、排列方式、坐标等)。.dat和.hdr文件共同构成了ENVI的标准栅格数据格式。当你在ENVI中打开一个.dat文件时，软件实际上是同时读取了.dat和.hdr文件来正确显示和处理图像。

（1）核心作用：存储像元值数据

.dat文件包含了图像最基本的元素——像元(像素)的数值本身。
这些数值代表了传感器捕获的地物辐射能量(通常经过一定处理，如辐射定标后的辐亮度、反射率等，或原始的DN值)。

（2）与.hdr文件的密不可分性

单独的.dat文件本身是"无意义"的，因为它不包含任何解释这些二进制数据含义的信息。
它必须与一个同名的.hdr(头文件) 文件配对使用。
.hdr文件是一个纯文本文件，包含了解读.dat文件所必需的元数据(Metadata)，例如：

samples：图像列数（宽度）
lines：图像行数（高度）。
bands：图像波段数。
data type：像元值的数据类型（如，1：无符号整数、2：有符号整数、4：浮点数等）。
interleave：数据的排列方式（BSQ-按波段顺序排列、BIL-按行交叉排列、BIP-按像元交叉排列）。
byte order：字节顺序（大端序1或小端序0）。
map info：地理坐标信息（投影、像素大小、左上角坐标等）。
wavelength：各个波段的中心波长（对于高光谱数据尤为重要）。
以及其他描述性信息（如传感器类型、采集时间等）。ENVI
description = {Band Math Result, Expression = [b1*0.0000275 - 0.2] B1:Band1:LC08_L2SP_124036_20240517_20240521_02_T1_SR_B2.TIF [Mon Apr 07 22:49:042025]}
samples = 7571
lines   = 7711
bands   = 1
header offset = 0
file type = ENVI Standard
data type = 4
interleave = bsq
sensor type = Unknown
byte order = 0
map info = {UTM, 1.000, 1.000, 609585.000, 3947715.000, 3.0000000000e+001, 3.0000000000e+001, 49, North, WGS-84, units=Meters}
coordinate system string = {PROJCS["WGS_1984_UTM_Zone_49N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",111.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]}
wavelength units = Unknown
band names = {Band Math (b1*0.0000275-0.2)}

（3）分离方式的优势

灵活性： 可以轻松修改元数据.hdr而不需要重新生成庞大的数据文件.dat。例如，可以快速更改投影信息或添加波长信息。
效率： 二进制存储.dat对于读写大型图像数据非常高效。
兼容性： 文本格式的.hdr易于人类阅读和编辑，也易于被不同软件解析(只要它们理解ENVI头文件格式)。
通用性： 这种"头文件+数据文件"的模式在科学数据处理领域非常常见。

（4）处理方式

可以参照"GeoTIFF"的处理方法，两者大同小异，基本一致。

⛄矢量数据

👀常用Python库

（1）GeoPandas

GeoPandas是用于地理数据处理和分析的Python库。它扩展了Pandas库，能够轻松处理地理数据。GeoPandas提供了丰富的功能，包括读取和写入地理数据文件、几何操作、空间连接、地理编码等。GeoPandas提供两个主要的数据结构：GeoSeries和GeoDataFrame。GeoSeries是一个扩展的Pandas Series对象，用于存储几何对象。GeoDataFrame是一个扩展的Pandas DataFrame对象，其中包含一个特殊的 geometry(GeoSeries)列，用于存储地理几何对象（如点、线、多边形等）。参考博客①

核心特点：

矢量数据操作(裁剪/合并/空间连接)
基于Pandas的数据结构(GeoSeries、GeoDataFrame)
坐标系管理(CRS转换)
空间可视化集成
支持多种文件格式(Shapefile/GeoJSON/GPKG)

主要依赖库：

Pandas(数据处理)
Shapely(几何操作)
Fiona(文件I/O)
Pyproj(坐标转换)

geopandas核心结构：

①GeoDataFrame是geopandas的核心数据结构，类似于pandas的DataFrame，但包含一个特殊的geometry列，用于存储地理几何对象。

②GeoSeries是geopandas的另一个核心数据结构，类似于pandas的Series，但专门用于存储地理几何对象。

geopandas常用函数：

①read_file()函数用于从文件中读取地理数据，支持多种格式，如Shapefile、GeoJSON、KML等。

②to_file()函数用于将GeoDataFrame或GeoSeries写入文件。

③plot()函数用于绘制地理数据。

④sjoin()函数用于空间连接(Spatial Join)，即将两个GeoDataFrame根据空间关系进行连接。连接方式(‘inner’，‘left’，‘right’)。

⑤overlay()函数用于执行空间叠加操作(Overlay)，如交集、并集、差异等。叠加方式(‘union’，‘difference’，‘intersection’，‘symmetric_difference’)。

⑥buffer()函数用于计算几何对象的缓冲区。

⑦centroid()函数用于计算几何对象的质心。

⑧dissolve()函数用于根据某一列的值对GeoDataFrame进行聚合。

⑨cx属性用于根据坐标范围筛选GeoDataFrame或GeoSeries。

…

（2）Shapely

Shapely是一个BSD许可的Python包，用于操作和分析笛卡尔平面中的几何对象。它使用广泛部署的开源几何库GEOS(PostGIS的引擎，JTS的一个端口)。Shapely封装了GEOS的几何形状和操作，为单个(标量)几何形状提供了丰富的几何接口，并为使用几何数组的操作提供了更高性能的NumPy ufuncs。

核心特点：

基础几何对象操作(点/线/面)
空间关系计算(相交/包含/距离)
几何变换(缓冲区/仿射变换)
拓扑操作(并集/交集/差分)

主要依赖库：

GEOS(C++几何引擎)
Numpy(数组支持)

（3）Fiona

Fiona是一个专为Python开发的地理空间矢量数据处理库，支持多种格式如Shapefile和GeoJSON。它通过简单的Python IO风格操作，帮助开发者读取和写入地理数据文件，同时与其他GIS工具(如GDAL、Shapely)无缝集成。Fiona的设计注重简洁和可靠，适合处理多层GIS格式数据。

核心特点：

多格式矢量数据读写(70+格式)
流式处理大文件
元数据访问
低级别要素操作

主要依赖库：

GDAL/OGR(地理数据抽象库)
Click(命令行支持)

（4）GDAL/OGR

OGR是GDAL的一个子项目，提供对矢量数据的支持。它实现了一个对空间参考信息进行处理的类，用来对空间数据的空间信息进行处理。GDAL提供了一整套读写不同栅格数据格式功能的抽象类库，而OGR是一个读写诸多矢量数据格式功能的抽象类库。两大类库用同一个生成系统进行维护，因此通常把这两个库合称为GDAL/OGR ，或者简称为GDAL 。

核心特点：

专业级空间数据转换
高级坐标系转换
格式互操作
栅格矢量互转

主要依赖库：

GDAL C库
Numpy

（5）Dask-geopandas

Dask-geopandas是一个结合了Dask和GeoPandas的框架，用于处理和分析大型地理空间数据。它通过并行计算和延迟执行来提高处理效率，特别适用于处理大规模的GIS数据。

核心特点：

分布式空间计算
大数据分块处理
延迟计算优化
集群部署支持

主要依赖库：

Dask(并行计算)
GeoPandas
Distributed(任务调度)

Dask库：Dask是一个开源的Python库，专为并行计算和大数据处理设计。它提供了与Pandas和NumPy类似的高层次接口，同时支持将计算分布到多核、集群或云环境中。Dask通过分块(chunking)和延迟计算(lazy evaluation)技术，实现了高效的数据处理和计算加速。核心功能：

并行数据处理：通过分块和延迟计算实现并行数据处理。
大数据处理：支持处理超出内存容量的大规模数据集。
数据帧操作：提供与Pandas类似的高层次数据帧接口。
数组操作：提供与NumPy类似的高层次数组接口。
并行计算：支持将计算任务分布到多核、集群或云环境中

👀Shapefile

（1）Shapefile矢量裁剪

import geopandas as gpdshp_path1 = r"D:\Desktop\data\shp\Layer1.shp"
shp_path2 = r"D:\Desktop\data\shp\Layer2.shp"
shp_clip_path = r"D:\Desktop\data\shp\clip_result.shp"# 读取被裁剪数据和裁剪范围数据
target = gpd.read_file(shp_path1)   # 被裁剪的矢量
clipper = gpd.read_file(shp_path2)  # 裁剪范围# 执行裁剪操作
clipped = gpd.clip(target, clipper)# 保存结果
clipped.to_file(shp_clip_path, encoding='utf-8')

（2）Shapefile矢量合并

# 简单叠加合并
import pandas as pd
import geopandas as gpd
import glob# 读取所有需要合并的Shapefile
files = glob.glob(r"D:\Desktop\data\shp\*.shp")
gdfs = [gpd.read_file(f) for f in files]# 确保所有几何类型一致（可选）
for gdf in gdfs:assert gdf.geom_type[0] in ('Polygon', 'MultiPolygon'), "几何类型不一致"# 合并数据
merged = gpd.GeoDataFrame(pd.concat(gdfs, ignore_index=True))# 保存结果
shp_merged_path = r"D:\Desktop\data\shp\merged_result.shp"
merged.to_file(shp_merged_path, encoding='utf-8')

# 拓扑合并（融合相邻边界）,只针对一个矢量图层
import geopandas as gpd# 读取数据
shp_merged_path = r"D:\Desktop\data\shp\merged_result.shp"
gdf = gpd.read_file(shp_merged_path)# 按属性字段融合边界（例如根据'province'字段合并）,一个图层中把多个面融合为一个面
# 将多行数据合并为一行，geometry合并，其他行只保留一行信息
# 如果需要按照多个列进行合并，可以传递一个列名的列表
dissolved = gdf.dissolve(by='Id', aggfunc='sum')
# dissolved = gdf.dissolve(by=['Id','Shape'], aggfunc='sum')# 无字段整体融合
unified = gdf.dissolve()# 保存结果
shp_dissolved_path = r"D:\Desktop\data\shp\dissolved_result.shp"
dissolved.to_file(shp_dissolved_path)

在Python中，assert语句用于测试一个表达式，如果表达式结果为False，则会触发AssertionError异常。这种方式可以在条件不满足程序运行的情况下直接返回错误，而不必等待程序运行后出现崩溃的情况。

（3）Shapefile矢量交集

# 矢量图层交集取反
import geopandas as gpdshp_path1 = r"D:\Desktop\data\shp\Layer1.shp"
shp_path2 = r"D:\Desktop\data\shp\Layer2.shp"gdf_left = gpd.read_file(shp_path1)
gdf_right = gpd.read_file(shp_path2)
# 计算矢量图层1和矢量图层2的交集部分，保留的geometry为矢量图层1去掉交集部分后的结果
gdf_left_diff_ritht = gpd.overlay(gdf_left, gdf_right,'difference')# 保存结果
shp_overlay_path = r"D:\Desktop\data\shp\overlay_difference_result.shp"
gdf_left_diff_ritht.to_file(shp_overlay_path)

# 矢量图层交集
import geopandas as gpdshp_path1 = r"D:\Desktop\data\shp\Layer1.shp"
shp_path2 = r"D:\Desktop\data\shp\Layer2.shp"gdf_left = gpd.read_file(shp_path1)
gdf_right = gpd.read_file(shp_path2)
# 计算矢量图层1和矢量图层2的交集部分，保留的geometry为交集部分
gdf_left_inte_ritht = gpd.overlay(gdf_left, gdf_right,'intersection')# 保存结果
shp_overlay_path = r"D:\Desktop\data\shp\overlay_intersection_result.shp"
gdf_left_inte_ritht.to_file(shp_overlay_path)

（4）区域统计

基于区域矢量数据，利用rasterstats库批量实现GTiff数据的区域统计，并将其值保存至CSV：

from rasterstats import zonal_stats
import pandas as pd
import os
import glob# 区域统计
def zonal_stats_simple(vector_path, raster_path, output_csv):raster_files = glob.glob(os.sep.join([raster_path, "*.tif"]))num = len(raster_files)result_df = pd.DataFrame()times = []mean_values = []for i in range(num):# 计算区域统计stats = zonal_stats(vector_path, raster_files[i], stats="mean", all_touched=False, nodata=0)time = os.path.basename(os.path.basename(raster_files[i]).split("_")[0])mean_value = [stat['mean'] for stat in stats]times.append(time)mean_values.append(mean_value[0])result_df["time"] = timesresult_df["mean_values"] = mean_values# 保存为CSVresult_df.to_csv(output_csv, index=False)if __name__ == "__main__":# 读取矢量数据vector_path = r"D:\Desktop\data\...\shp\Sample.shp"raster_path = r"D:\Desktop\data\...\2024_ERA5\rh1000hpa"output_csv = r"D:\Desktop\data\...\excel\rh1000hpa.csv"zonal_stats_simple(vector_path, raster_path, output_csv)

‌rasterstats是一个用于处理栅格和矢量数据的Python库，主要用于从栅格数据中提取基于矢量区划的统计信息，包括区域统计、分类统计和面状统计等，简化了栅格与矢量数据的交互‌。主要功能包括：

‌zonal_stats：基于矢量区域从栅格数据中提取统计信息，如均值、最小值、最大值、方差等‌。

‌point_query：用于从栅格中提取给定点的数值信息‌。

gen_zonal_stats：生成基于多边形的统计结果，适用于大规模的并行处理‌

👀GeoJSON

GeoJSON是一种对各种地理数据结构进行编码的格式，基于Javascript对象表示法(JSON)的地理空间信息数据交换格式。GeoJSON对象可以表示几何、特征或者特征集合。GeoJSON支持的几何类型包括点、线、面、多点、多线、多面和几何集合。GeoJSON里的特征包含一个几何对象和其他属性，特征集合表示一系列特征。

{"type": "FeatureCollection","features": [{"type": "Feature","properties": {},"geometry": {"coordinates": [[[113.68846316466124, 34.83541413400329], [113.68846316466124, 34.79070319327178], [113.76913047352139, 34.79070319327178], [113.76913047352139, 34.83541413400329], [113.68846316466124, 34.83541413400329]]],"type": "Polygon"}}]
}

利用GeoPandas库处理GeoJSON格式数据，并将其保存为Shapefile格式的代码示例：

import geopandas as gpd# 1. 读取GeoJSON文件（假设为面数据）
input_geojson = r"D:\Desktop\data\geojson\longhu.geojson"
gdf = gpd.read_file(input_geojson)# 检查数据是否为面类型（可选）
if not all(gdf.geometry.type.isin(['Polygon', 'MultiPolygon'])):print("警告：数据中包含非面类型几何体！")# GeoPandas会保留原始数据的坐标参考系统（CRS）。如果需要强制转换，可添加：
# gdf = gdf.to_crs("EPSG:4326")# 2. 保存为Shapefile
# Shapefile实际由多个文件组成（.shp、.shx、.dbf等），GeoPandas 会自动生成。
output_shapefile = r"D:\Desktop\data\geojson\output_data.shp"
gdf.to_file(output_shapefile, driver="ESRI Shapefile")

👀KML/KMZ

KML(Keyhole Markup Language)是一种基于XML的地理数据标记语言，由Google旗下的Keyhole公司开发(后成为Google Earth的核心格式)，主要用于存储点、线、面等地理要素的存储和可视化，尤其适合用于三维轨迹数据展示和地图标注。‌‌KML的核心特点：

XML结构：KML文件本质上是XML文件，结构清晰，标签语义明确，易于阅读和编写。‌‌
三维可视化：支持三维模型和效果(如拉伸多边形)，适合展示复杂的地理空间数据。‌‌
跨平台兼容性：KML是开放标准，支持多种GIS平台(如Google Earth、ArcGIS、QGIS等)，可在不同工具间轻松导入和导出。‌‌
轻量级：文件通常较小，适合网络传输和在线加载。

‌KML的应用场景：

地图标注与路径规划：如旅游路线标注、无人机航线生成等。‌‌
三维数据展示：适用于毕业设计或项目中展示三维立体图。‌‌
数据共享：可通过KMZ(压缩版KML)打包分发包含图片、图标等资源的复杂地理数据。‌‌

KMZ(Keyhole Markup Zip)：KML的压缩版本，本质是一个ZIP压缩包(扩展名为.kmz)，其用途是将KML文件及其关联资源(如图片、模型、图标等)打包为单一文件。其特点：

节省存储空间，便于共享。
必须包含至少一个doc.kml(主文件)，其他资源(如纹理、Collada模型)存储在压缩包内。

特性	KML	KMZ
格式	`XML`文本文件	`ZIP`压缩包(内含`KML`+资源)
扩展名	`.KML`	`.KMZ`
文件大小	较大(纯文本)	较小(压缩后)
资源支持	需外部链接图片等	可内嵌图片、模型等
编辑性	可直接编辑	需解压后编辑

<kml xmlns="http://www.opengis.net/kml/2.2"><Document><Placemark><ExtendedData/>      <!-- 空标签 --><Polygon><outerBoundaryIs><LinearRing><coordinates>113.68846316466124,34.83541413400329113.68846316466124,34.79070319327178113.76913047352139,34.79070319327178113.76913047352139,34.83541413400329113.68846316466124,34.83541413400329</coordinates></LinearRing></outerBoundaryIs></Polygon></Placemark></Document>
</kml>

利用GeoPandas库处理KML格式数据，并将其保存为Shapefile格式的代码示例：

import geopandas as gpd
import fiona# 利用fiona库，激活KML文件读取格式
fiona.drvsupport.supported_drivers['KML'] = 'rw'# 读取KML文件
input_kml = r"D:\Desktop\data\kml\longhu.kml"
gdf = gpd.read_file(input_kml, driver='KML')
print(gdf.head())# 保存为Shapefile
output_shapefile = r"D:\Desktop\data\kml\output_data.shp"
gdf.to_file(output_shapefile, driver="ESRI Shapefile")

（1）KMZ文件：需先解压提取内部的KML，利用geopandas读取解压后的KML(通常为doc.kml)。

（2）pykml库解析复杂KML：若需提取样式、描述等元数据，可使用pykml库。

kml文件中<ExtendedData>是一个非常重要的标签，用于存储自定义属性或元数据，这些数据可以附加到KML 的要素(如Placemark、Polygon等)上。它类似于Shapefile的"属性表"或GeoJSON的properties字段，但提供了更灵活的存储方式。<ExtendedData>标签的作用：

存储额外的属性数据
- 例如：地名、描述、统计信息、ID、时间戳等非几何信息。
兼容Google Earth和其他GIS工具
- Google Earth会解析这些数据并在“详细信息”面板中显示。
支持多种数据格式
- 可以存储文本、数字、嵌套结构(如JSON)，甚至二进制数据(需编码)。

<!-- 示例 -->
<Placemark><name>北京天安门</name><ExtendedData><Data name="population"><displayName>人口</displayName><value>21540000</value></Data><Data name="founded"><displayName>建城时间</displayName><value>1420</value></Data></ExtendedData><Point><coordinates>116.3975,39.9087</coordinates></Point>
</Placemark>

👀WKT

‌WKT(Well-Known Text)是一种用于描述地理空间几何对象的文本格式标准，由开放地理空间联盟OGC定义和维护。WKT是一种基于文本的地理空间数据表示方法，通过字符序列描述点、线、多边形等几何对象的形状和空间关系。‌‌常见几何类型包括点(Point)、线(LineString)、面(Polygon)、多点(MultiPoint)、复合面(MultiPolygon)。

简洁易读，适合存储单个几何对象。
不支持属性数据(仅描述几何形状)。
扩展格式：EWKT(PostGIS扩展，支持SRID信息)。

WKT示例数据如下(第一行为面数据，第二行为点数据)：

POLYGON ((113.68846316466124 34.83541413400329, 113.68846316466124 34.79070319327178, 113.76913047352139 34.79070319327178, 113.76913047352139 34.83541413400329, 113.68846316466124 34.83541413400329))
POINT (113.72885639495456 34.8140740745182)

（1）利用GeoPandas和shapely库处理WKT格式数据，并将其保存为Shapefile格式的代码示例一：

import geopandas as gpd
from shapely.wkt import loads# 1. 读取WKT文件（假设每行一个WKT面数据）
wkt_file = r"D:\Desktop\data\wkt\longhu_polygon_point.wkt"
output_shp = "D:\Desktop\data\wkt\output_polygons.shp"
# 2. 读取并解析WKT数据（每行一个WKT字符串）
with open(wkt_file, 'r') as f:wkt_lines = [line.strip() for line in f if line.strip()]
print(wkt_lines)
# 3. 创建几何对象列表和属性表
geometries = []
properties = []
for i, wkt in enumerate(wkt_lines, start=1):try:geom = loads(wkt)# 确保是面数据（Polygon或MultiPolygon）if geom.type in ('Polygon', 'MultiPolygon'):geometries.append(geom)properties.append({'id': i, 'wkt': wkt[:]})  # 示例属性else:print(f"警告：跳过非面几何体（第{i}行）：{geom.type}")except Exception as e:print(f"解析错误（第{i}行）：{e}")
# 4. 创建GeoDataFrame
if geometries:gdf = gpd.GeoDataFrame(properties, geometry=geometries, crs="EPSG:4326")# 5. 保存为Shapefilegdf.to_file(output_shp, driver="ESRI Shapefile")print(f"成功保存Shapefile：{output_shp}")print(f"包含{len(gdf)}个面要素")
else:print("错误：未找到有效的面数据")# 输出结果
['POLYGON ((113.68846316466124 34.83541413400329, 113.68846316466124 34.79070319327178, 113.76913047352139 34.79070319327178, 113.76913047352139 34.83541413400329, 113.68846316466124 34.83541413400329))', 'POINT (113.72885639495456 34.8140740745182)']
警告：跳过非面几何体（第2行）：Point
成功保存Shapefile：D:\Desktop\data\wkt\output_polygons.shp
包含1个面要素

（2）利用GeoPandas和shapely库处理WKT格式数据，并将其保存为Shapefile格式的代码示例二：

import geopandas as gpd
from shapely.wkt import loads# 示例WKT面数据（Polygon 或 MultiPolygon）
wkt_data = ["POLYGON ((116.404 39.915, 116.404 39.920, 116.410 39.920, 116.410 39.915, 116.404 39.915))","POLYGON ((116.415 39.925, 116.415 39.930, 116.420 39.930, 116.420 39.925, 116.415 39.925))"
]
# 示例属性数据（可选）
attributes = [{"id": 1, "name": "区域A"},{"id": 2, "name": "区域B"}
]
# 1. 将WKT字符串转换为几何对象
geometries = [loads(wkt) for wkt in wkt_data]
# 2. 创建GeoDataFrame
gdf = gpd.GeoDataFrame(attributes, geometry=geometries, crs="EPSG:4326")  # 默认WGS84坐标系
# 3. 保存为Shapefile
output_path = r"D:\Desktop\data\wkt\output_polygons2.shp"
gdf.to_file(output_path, driver="ESRI Shapefile", encoding="UTF-8")
print(f"Shapefile 已保存至：{output_path}")# 输出结果
Shapefile 已保存至：D:\Desktop\data\wkt\output_polygons2.shp

👀GML

‌GML(Geography Markup Language)是一种基于XML的地理信息编码语言，由开放地理空间联盟OGC制定，主要用于地理数据的建模、存储和交换，旨在解决不同来源、模型和格式的地理数据共享与互操作问题。其核心功能包括：‌‌

编码地理要素的几何对象和属性信息，支持复杂空间要素和拓扑关系。
实现内容与表现的分离，支持灵活的数据解析与可视化。

GML示例数据如下:

<gml:FeatureCollection xmlns:gml="http://www.opengis.net/gml/3.2"><gml:featureMember><gml:Polygon gml:id="polygon1"><gml:extentOf><gml:Polygon srsName="urn:ogc:def:crs:epsg::4326"><gml:exterior><gml:LinearRing><gml:posList>40.071713788991 116.374131643854 40.0715518439545 116.375167866672 40.0710444161736 116.375735659996 40.0703318580131 116.375636296165 40.0699863752687 116.37506850284 40.0700727459548 116.373691604027 40.0708824711371 116.373052836537 40.0709688418232 116.37302444687 40.0713791025823 116.3732515642 40.071713788991 116.374131643854</gml:posList></gml:LinearRing></gml:exterior></gml:Polygon></gml:extentOf></gml:Polygon></gml:featureMember>
</gml:FeatureCollection>

（1）利用GeoPandas库处理GML格式数据，并将其保存为GeoJSON格式的代码示例一：

import geopandas as gpd# 直接读取GML文件
input_gml = r"D:\Desktop\data\gml\GML_sample.gml"
output_geojson = r"D:\Desktop\data\gml\output.geojson"
gdf = gpd.read_file(input_gml, driver="GML")
gdf.to_file(output_geojson, driver="GeoJSON")

（2）利用GDAL库的ogr2ogr处理GML格式数据，并将其保存为Shapefile格式的代码示例二：

import subprocessdef readGML(gml_path, shp_path):# 原始GML文件的投影信息的格式可能存在问题，不能匹配cmd = ["ogr2ogr","-f", "ESRI Shapefile","-t_srs", "EPSG:4326",shp_path,gml_path]subprocess.run(cmd, check=True)if __name__ == "__main__":input_gml = r"D:\Desktop\data\gml\GML_sample.gml"output_shp = r"D:\Desktop\data\gml\output.shp"readGML(input_gml, output_shp)

subprocess模块：subprocess模块是Python标准库中的一个模块，用于创建和管理子进程，允许用户执行外部命令、连接到输入/输出流，并获取子进程的返回码。该模块的主要功能包括执行外部命令、控制输入/输出、错误处理以及进程管理。

ogr2ogr工具：‌ ogr2ogr是一个开源的命令行工具，属于GDAL库的一部分，主要用于处理和转换地理空间数据。它支持多种矢量数据格式的转换，包括Shapefile、GeoJSON、KML、GML、GeoPackage等，是地理信息系统数据处理中不可或缺的工具‌。

⛄气象数据格式

👀NetCDF

NetCDF(network Common Data Form)网络通用数据格式是由美国大学大气研究协会(UCAR)的Unidata项目科学家针对科学数据的特点开发的，是一种面向数组型并适于网络共享的数据的描述和编码标准。NC格式是一种用于存储多维科学数据的通用文件格式，广泛应用于气候科学、地球科学、水文等领域。它支持数组型数据，包含变量(Variables)、维度(Dimensions)和属性(Attributes)三种核心结构，便于网络共享和跨平台处理。NetCDF文件的数据结构包含三个主要组成部分：

‌变量（Variables）：存储实际的多维数组数据(如温度、降水等)，每个变量关联到特定的维度。
‌维度（Dimensions）：定义变量的维度信息(如时间、经度、纬度等)，类似于变量的坐标轴。‌‌
‌属性（Attributes）：存储元数据信息。

利用netCDF4和GDAL库处理NetCDF格式数据，并将其保存为GTiff格式的代码示例：

格式转换：NetCDF→GTiff
逐小时数据自动进行时区转换（可选）

import netCDF4 as nc
from osgeo import gdal, osr
import numpy as np
import os
import cftime
import pandas as pd
from pytz import timezone# 写入影像
def write_img(filename, im_proj, im_geotrans, im_data):# 判断栅格数据的数据类型if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# 判读数组维数if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelse:im_bands, (im_height, im_width) = 1, im_data.shape# 创建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数dataset.SetProjection(im_proj) # 写入投影if im_bands == 1:dataset.GetRasterBand(1).WriteArray(im_data) # 写入数组数据else:for i in range(im_bands):dataset.GetRasterBand(i + 1).WriteArray(im_data[i,:,:])del dataset# 将NetCDF数据转化为GTiff
def nc_totif(input_path, output_path):# 读取nc文件tep_data = nc.Dataset(input_path)print(tep_data.variables.keys())# 获取nc文件中对应变量的信息lon_data = tep_data.variables["longitude"][:]lat_data = tep_data.variables["latitude"][:]# 影像的左上角&右下角坐标lonmin, latmax, lonmax, latmin = [lon_data.min(), lat_data.max(), lon_data.max(), lat_data.min()]# 分辨率计算num_lon = len(lon_data)num_lat = len(lat_data)lon_res = (lonmax - lonmin) / (float(num_lon) - 1)lat_res = (latmax - latmin) / (float(num_lat) - 1)# 定义投影proj = osr.SpatialReference()proj.ImportFromEPSG(4326) 		# WGS84proj = proj.ExportToWkt() 		# 重点，转成wkt格式# 仿射变换矩阵geotransform = (lonmin, lon_res, 0.0, latmax, 0.0, -lat_res)# 获取时间time_var = tep_data.variables["valid_time"]time_units = time_var.unitsprint(time_units)time_values = time_var[:]# 将时间戳转换为日期时间对象datetimes = cftime.num2date(time_values, units=time_units, calendar="gregorian")dates = [dt.strftime("%Y%m%d") for dt in datetimes]times = [dt.strftime("%H%M%S") for dt in datetimes]dates = np.unique(dates)'''# 获取时间(逐小时一般进行时区转换)time_var = tep_data.variables["valid_time"]  # 获取时间变量time_units = time_var.units   # 获取时间单位print(time_units)time_values = time_var[:]     # 获取时间值# 将时间戳转换为日期时间对象datetimes = cftime.num2date(time_values, units = time_units, calendar = "gregorian")utc_time_str = [dt.strftime("%Y-%m-%d %H:%M:%S") for dt in datetimes]beijing_time_str = []for i in range(len(utc_time_str)):utc_time = pd.to_datetime(utc_time_str[i], format="%Y-%m-%d %H:%M:%S")utc_time = utc_time.tz_localize("UTC")beijing_time = utc_time.astimezone(timezone("Asia/Shanghai")).strftime("%Y%m%d%H")beijing_time_str.append(beijing_time)beijing_time_str = np.unique(beijing_time_str)print(beijing_time_str)'''t2m = tep_data.variables["t2m"][:]t2m_arr = np.asarray(t2m)t2m_day_mean = np.mean(t2m_arr, 0)outpath = os.sep.join([output_path, dates[0]+"_daymean_t2m.tif"])write_img(outpath, proj, geotransform, t2m_day_mean)if __name__ == "__main__":# nc文件输入输出路径input_path = r"D:\Desktop\data\20250528\20250528_ERA5_t2m_wind_tp.nc"output_path = r"D:\Desktop\data\20250528\tif"# 读取nc文件，转换为tif文件nc_totif(input_path, output_path)

👀GRIB

GRIB(Gridded Binary)是世界气象组织(WMO)开发的一种用于存储和传输网格化气象数据的二进制文件格式，广泛应用于数值天气预报、气候研究等领域。其基本特点：

二进制格式：紧凑高效，适合大量数据存储和传输
自描述性：包含元数据描述数据内容
压缩存储：支持多种压缩算法
分节结构：数据按节(Section)组织
国际标准：WMO标准格式，全球气象业务通用

import xarray as xr
import rioxarray
from pathlib import Path
import warningsdef grib_to_tiff(grib_path, tiff_path, variable=None, time_idx=0, level_idx=0):"""使用xarray和cfgrib将GRIB文件转换为GeoTIFF参数:grib_path (str): 输入GRIB文件路径tiff_path (str): 输出TIFF文件路径variable (str): 可选，指定要提取的变量名time_idx (int): 可选，时间维度索引(默认为0)level_idx (int): 可选，高度层索引(默认为0)"""try:# 禁用不必要的警告warnings.filterwarnings("ignore", category=UserWarning)# 打开GRIB文件ds = xr.open_dataset(grib_path, engine="cfgrib")# 如果没有指定变量，选择第一个数据变量if variable is None:variable = list(ds.data_vars.keys())[0]print(f"未指定变量，自动选择: {variable}")# 获取数据数组data = ds[variable]# 处理多维数据(时间、高度层等)if "time" in data.dims:data = data.isel(time=time_idx)if "isobaricInhPa" in data.dims:data = data.isel(isobaricInhPa=level_idx)if "heightAboveGround" in data.dims:data = data.isel(heightAboveGround=level_idx)# 添加地理参考信息data = data.rio.set_spatial_dims(x="longitude", y="latitude", inplace=False)data.rio.write_crs("EPSG:4326", inplace=True)  # WGS84坐标系# 保存为GeoTIFFdata.rio.to_raster(tiff_path)print(f"成功转换: {Path(grib_path).name} → {Path(tiff_path).name}")return Trueexcept Exception as e:print(f"转换失败: {str(e)}")return Falsefinally:warnings.resetwarnings()# 使用示例
if __name__ == "__main__":input_file = "example.grib"  # 输入GRIB文件output_file = "output.tif"   # 输出TIFF文件# 基本转换grib_to_tiff(input_file, output_file)# 高级用法：指定变量和时间层# 850hPa温度# grib_to_tiff(input_file, "temp_850hPa.tif", variable="t", level_idx=2)

GRIB数据读取方法：

（1）使用xarray+cfgrib组合

（2）使用pygrib库，小编一直安装错误…

NetCDF和GRIB是气象和地球科学领域常用的两种数据格式，它们的主要区别如下：

特性	NetCDF	GRIB
设计目的	通用科学数据格式	气象领域专用格式
标准组织	UCAR/Unidata	WMO(世界气象组织)
版本	NetCDF3/4	GRIB1/GRIB2
数据结构	多维数组+元数据	消息集合(每个消息独立)
压缩效率	中等(支持压缩)	高(专为气象数据优化)
元数据灵活性	高(可自定义属性)	固定模板(参数代码表)
时间处理	灵活的时间坐标	基于参考时间+预报时长
空间参考	可自由定义	预定义投影系统
编程接口	丰富(多种语言支持)	较复杂(需专用库)
典型扩展名	`.nc`、`.cdf`	`.grb`、`.grib`
主要用途	科研数据存储/交换	气象业务数据传输