栅格数据以规则网格(像素)的数值矩阵表达地理现象，每个单元格代表一个属性值(如高程、温度)。例如卫星影像、数字高程模型、温度分布图。存储格式包括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(核心依赖)NumPyaffine(变换矩阵)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库NumPympi4py(可选，并行支持)
👀GeoTIFF数据
（1）GDAL方式
影像读取
影像写入
基于矢量边界的影像裁剪Warp()。参考博客① 、参考博客②
虚拟内存读写数据（文件路径前缀添加/vsimem/）
......
from osgeo import gdaldef read_img (filename ): dataset = gdal.Open(filename) im_Description = dataset.GetDescription() im_band = dataset.GetRasterBand(1 ) print (im_band.GetMinimum()) print (im_band.GetMaximum()) im_width = dataset.RasterXSize im_height = dataset.RasterYSize im_bands = dataset.RasterCount im_geotrans = dataset.GetGeoTransform() im_proj = dataset.GetProjection() im_data = dataset.ReadAsArray(0 , 0 , im_width, im_height).astype(float ) del dataset return im_proj, im_geotrans, im_data, im_height, im_width
from osgeo import gdaldef write_img (filename, im_proj, im_geotrans, im_data ): if 'int8' in im_data.dtype.name: datatype = gdal.GDT_Byte elif 'int16' in im_data.dtype.name: datatype = gdal.GDT_UInt16 else : datatype = gdal.GDT_Float32 if len (im_data.shape) == 3 : im_bands, im_height, im_width = im_data.shape else : im_bands, (im_height, im_width) = 1 , im_data.shape driver = gdal.GetDriverByName("GTiff" ) options = ['INTERLEAVE=BAND' ] 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 Pixel dataset = 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 osimport rasteriotemp_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) 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 rasteriofilename = 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 rasterioimport geopandas as gpdfrom rasterio.mask import maskdef 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 globimport osimport rasteriofrom rasterio.merge import mergedef 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" 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 gpdimport rasteriofrom rasterio.io import MemoryFilefrom rasterio.mask import maskdef 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 )) 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: 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_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的区别：
返回类型 字典(可写，可直接用于创建新文件)
字典(只读副本)
主要用途 创建新栅格时的模板
查看元数据
可变性 可修改(修改后会影响数据集)
不可修改
包含内容 核心元数据 + 部分驱动特定选项
仅核心元数据
版本引入 Rasterio 1.0+早期版本即有
rasterio中dataset.meta属性设置及其详细说明：
driverstr栅格驱动名称(如'GTiff'、'PNG'、'HFA'等)
dtypestr数据类型(如'uint8'、'float32'等)，对应NumPy数据类型
nodataint/float/None表示无效数据的值(如0、-9999或None表示无Nodata值)
widthint栅格列数(宽度)
heightint栅格行数(高度)
countint波段数量(如RGB图像为3)
crsrasterio.crs.CRS坐标参考系统(如CRS.from_epsg(4326)表示WGS84)
transformAffine仿射变换矩阵
{ '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 ) } { '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。
import osimport h5pyimport numpy as npfrom osgeo import gdal, osrdef write_img (filename, im_proj, im_geotrans, im_data ): if 'int8' in im_data.dtype.name: datatype = gdal.GDT_Byte elif 'int16' in im_data.dtype.name: datatype = gdal.GDT_UInt16 else : datatype = gdal.GDT_Float32 if len (im_data.shape) == 3 : im_bands, im_height, im_width = im_data.shape else : 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 def hdf_to_tif (input_folder_path, output_path ): 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) ''' HDF5_group = f["HDFEOS" ]["GRIDS" ]["VIIRS_Grid_DNB_2d" ]["Data Fields" ] DNB_data = HDF5_group["AllAngle_Composite_Snow_Free" ][:]*0.1 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 ) proj = proj.ExportToWkt() geotransform = (lonmin, lon_res, 0.0 , latmax, 0.0 , -lat_res) 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__' : input_folder_path = r"D:\Desktop\data\hdf\VNP46A4" output_path = r"D:\Desktop\data\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:Band 1:LC08_L2SP_124036_20240517_20240521_02_T1_SR_B2.TIF [Mon Apr 07 22:49:04 2025]} 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 。
核心特点 ：
专业级空间数据转换
高级坐标系转换
格式互操作
栅格矢量互转
主要依赖库 ：
（5）Dask-geopandas
Dask-geopandas是一个结合了Dask和GeoPandas的框架，用于处理和分析大型地理空间数据。它通过并行计算和延迟执行来提高处理效率，特别适用于处理大规模的GIS数据。
核心特点 ：
分布式空间计算
大数据分块处理
延迟计算优化
集群部署支持
主要依赖库 ：
Dask(并行计算)GeoPandasDistributed(任务调度)
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 pdimport geopandas as gpdimport globfiles = 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 gpdshp_merged_path = r"D:\Desktop\data\shp\merged_result.shp" gdf = gpd.read_file(shp_merged_path) dissolved = gdf.dissolve(by='Id' , 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) 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) 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_statsimport pandas as pdimport osimport globdef 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" ] = times result_df["mean_values" ] = mean_values result_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 gpdinput_geojson = r"D:\Desktop\data\geojson\longhu.geojson" gdf = gpd.read_file(input_geojson) if not all (gdf.geometry.type .isin(['Polygon' , 'MultiPolygon' ])): print ("警告：数据中包含非面类型几何体！" ) 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模型)存储在压缩包内。
格式 XML文本文件ZIP压缩包(内含KML+资源)
扩展名 .KML.KMZ
文件大小 较大(纯文本)
较小(压缩后)
资源支持 需外部链接图片等
可内嵌图片、模型等
编辑性 可直接编辑
需解压后编辑
<kml xmlns ="http://www.opengis.net/kml/2.2" > <Document > <Placemark > <ExtendedData /> <Polygon > <outerBoundaryIs > <LinearRing > <coordinates > 113.68846316466124,34.83541413400329 113.68846316466124,34.79070319327178 113.76913047352139,34.79070319327178 113.76913047352139,34.83541413400329 113.68846316466124,34.83541413400329 </coordinates > </LinearRing > </outerBoundaryIs > </Polygon > </Placemark > </Document > </kml >
利用GeoPandas库处理KML格式数据，并将其保存为Shapefile格式的代码示例：
import geopandas as gpdimport fionafiona.drvsupport.supported_drivers['KML' ] = 'rw' input_kml = r"D:\Desktop\data\kml\longhu.kml" gdf = gpd.read_file(input_kml, driver='KML' ) print (gdf.head())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 gpdfrom shapely.wkt import loadswkt_file = r"D:\Desktop\data\wkt\longhu_polygon_point.wkt" output_shp = "D:\Desktop\data\wkt\output_polygons.shp" with open (wkt_file, 'r' ) as f: wkt_lines = [line.strip() for line in f if line.strip()] print (wkt_lines)geometries = [] properties = [] for i, wkt in enumerate (wkt_lines, start=1 ): try : geom = loads(wkt) 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} " ) if geometries: gdf = gpd.GeoDataFrame(properties, geometry=geometries, crs="EPSG:4326" ) gdf.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 gpdfrom shapely.wkt import loadswkt_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" } ] geometries = [loads(wkt) for wkt in wkt_data] gdf = gpd.GeoDataFrame(attributes, geometry=geometries, crs="EPSG:4326" ) 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 gpdinput_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 ): 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 ncfrom osgeo import gdal, osrimport numpy as npimport osimport cftimeimport pandas as pdfrom pytz import timezonedef write_img (filename, im_proj, im_geotrans, im_data ): if 'int8' in im_data.dtype.name: datatype = gdal.GDT_Byte elif 'int16' in im_data.dtype.name: datatype = gdal.GDT_UInt16 else : datatype = gdal.GDT_Float32 if len (im_data.shape) == 3 : im_bands, im_height, im_width = im_data.shape else : 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 def nc_totif (input_path, output_path ): tep_data = nc.Dataset(input_path) print (tep_data.variables.keys()) 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 ) proj = proj.ExportToWkt() geotransform = (lonmin, lon_res, 0.0 , latmax, 0.0 , -lat_res) 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" ) 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__" : input_path = r"D:\Desktop\data\20250528\20250528_ERA5_t2m_wind_tp.nc" output_path = r"D:\Desktop\data\20250528\tif" nc_totif(input_path, output_path)
👀GRIB
GRIB(Gridded Binary)是世界气象组织(WMO)开发的一种用于存储和传输网格化气象数据的二进制文件格式，广泛应用于数值天气预报、气候研究等领域。其基本特点：
二进制格式 ：紧凑高效，适合大量数据存储和传输自描述性 ：包含元数据描述数据内容压缩存储 ：支持多种压缩算法分节结构 ：数据按节(Section)组织国际标准 ：WMO标准格式，全球气象业务通用
import xarray as xrimport rioxarrayfrom pathlib import Pathimport 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) 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 ) data.rio.to_raster(tiff_path) print (f"成功转换: {Path(grib_path).name} → {Path(tiff_path).name} " ) return True except Exception as e: print (f"转换失败: {str (e)} " ) return False finally : warnings.resetwarnings() if __name__ == "__main__" : input_file = "example.grib" output_file = "output.tif" grib_to_tiff(input_file, output_file)
GRIB数据读取方法：
（1）使用xarray+cfgrib组合
（2）使用pygrib库，小编一直安装错误......
NetCDF和GRIB是气象和地球科学领域常用的两种数据格式，它们的主要区别如下：
设计目的 通用科学数据格式
气象领域专用格式
标准组织 UCAR/Unidata
WMO(世界气象组织)
版本 NetCDF3/4
GRIB1/GRIB2
数据结构 多维数组+元数据
消息集合(每个消息独立)
压缩效率 中等(支持压缩)
高(专为气象数据优化)
元数据灵活性 高(可自定义属性)
固定模板(参数代码表)
时间处理 灵活的时间坐标
基于参考时间+预报时长
空间参考 可自由定义
预定义投影系统
编程接口 丰富(多种语言支持)
较复杂(需专用库)
典型扩展名 .nc、.cdf.grb、.grib
主要用途 科研数据存储/交换
气象业务数据传输
wechat(微信)
