基于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`
主要用途	科研數據存儲/交換	氣象業務數據傳輸