python 按像素提取固定经纬度的矩形区域栅格tif数据的值
我有栅格数据是已经处理好的数据，经过了几何校正辐射定标啥的
# -*- coding:utf-8 -*- """ @author: suyue @file: lianxi.py @time: 2025/09/09 @desc: 提取栅格数据的值 """ import rasterio from rasterio.mask import mask import numpy as np import pandas as pd import shapely from shapely.geometry import box import json def extract_raster_to_excel(raster_path, output_excel_path, min_lon, max_lon, min_lat, max_lat, resolution_meters=4000): """ 提取指定经纬度范围内的栅格数据并导出到Excel 参数: raster_path: 输入栅格文件路径 output_excel_path: 输出Excel文件路径 min_lon, max_lon: 经度范围 min_lat, max_lat: 纬度范围 resolution_meters: 栅格分辨率（米），用于估算网格 """ try: # 打开栅格文件 with rasterio.open(raster_path) as src: print(f"成功打开栅格文件: {raster_path}") print(f"栅格形状: {src.shape}") print(f"栅格范围: {src.bounds}") print(f"坐标参考系统: {src.crs}") # 创建矩形几何体 geom = box(min_lon, min_lat, max_lon, max_lat) geom_dict = [json.loads(shapely.to_geojson(geom))] # 提取指定范围内的数据 out_image, out_transform = mask(src, geom_dict, crop=True, all_touched=True) # 获取无数据值 nodata = src.nodata # 将三维数组转换为二维（如果是单波段） if out_image.ndim == 3: out_image = out_image[0] # 创建坐标网格 height, width = out_image.shape rows, cols = np.meshgrid(np.arange(height), np.arange(width), indexing='ij') # 计算每个像元的中心点坐标 lon_coords = [] lat_coords = [] for row in range(height): for col in range(width): # 计算像元中心点的坐标 lon, lat = rasterio.transform.xy(out_transform, row + 0.5, col + 0.5) lon_coords.append(lon) lat_coords.append(lat) # 准备数据 data = { '行号': rows.flatten(), '列号': cols.flatten(), '经度': lon_coords, '纬度': lat_coords, '栅格值': out_image.flatten() } # 创建DataFrame df = pd.DataFrame(data) # 过滤掉无数据值 if nodata is not None: df = df[df['栅格值'] != nodata] # 保存到Excel df.to_excel(output_excel_path, index=False, engine='openpyxl') print(f"成功提取数据到: {output_excel_path}") print(f"提取了 {len(df)} 个有效像元") print(f"数据范围: 经度 {df['经度'].min():.6f} - {df['经度'].max():.6f}") print(f" 纬度 {df['纬度'].min():.6f} - {df['纬度'].max():.6f}") return df except Exception as e: print(f"处理过程中出现错误: {e}") return None # 使用示例 if __name__ == "__main__": # 您的参数设置 input_raster = "D:/20240809example/CTT/CTT_20240810_0700.tif" # 替换为您的栅格文件路径 output_excel = "D:/20240809example/CTT/CTT_20240810_0700提取结果.xlsx" # 您指定的经纬度范围 min_longitude = 114.2692 max_longitude = 117.2414 min_latitude = 43.2313 max_latitude = 45.32772 # 执行提取 result_df = extract_raster_to_excel( raster_path=input_raster, output_excel_path=output_excel, min_lon=min_longitude, max_lon=max_longitude, min_lat=min_latitude, max_lat=max_latitude, resolution_meters=4000 ) # 显示前几行数据 if result_df is not None: print("\n前5行数据预览:") print(result_df.head())
输出
成功打开栅格文件: D:/20240809example/CTT/CTT_20240810_0700.tif
栅格形状: (612, 1250)
栅格范围: BoundingBox(left=85.0, bottom=31.968, right=130.0, top=54.0)
坐标参考系统: EPSG:4326
成功提取数据到: D:/20240809example/CTT/CTT_20240810_0700提取结果.xlsx
提取了 4980 个有效像元
数据范围: 经度 114.304000 - 117.256000
纬度 43.200000 - 45.324000
前5行数据预览:
行号 列号 经度 纬度 栅格值
0 0 0 114.304 45.324 -27.757851
1 0 1 114.340 45.324 -27.155380
2 0 2 114.376 45.324 -26.854853
3 0 3 114.412 45.324 -27.300863
4 0 4 114.448 45.324 -25.919701
Process finished with exit code 0
批量处理数据代码
# -*- coding:utf-8 -*- """ @author: suyue @file: lianxi.py @time: 2025/09/09 @desc: 批量提取栅格数据的值 """ import rasterio from rasterio.mask import mask import numpy as np import pandas as pd import shapely from shapely.geometry import box import json import os import glob def extract_raster_to_excel(raster_path, output_excel_path, min_lon, max_lon, min_lat, max_lat, resolution_meters=4000): """ 提取指定经纬度范围内的栅格数据并导出到Excel 参数: raster_path: 输入栅格文件路径 output_excel_path: 输出Excel文件路径 min_lon, max_lon: 经度范围 min_lat, max_lat: 纬度范围 resolution_meters: 栅格分辨率（米），用于估算网格 """ try: # 打开栅格文件 with rasterio.open(raster_path) as src: print(f"成功打开栅格文件: {raster_path}") print(f"栅格形状: {src.shape}") print(f"栅格范围: {src.bounds}") print(f"坐标参考系统: {src.crs}") # 创建矩形几何体 geom = box(min_lon, min_lat, max_lon, max_lat) geom_dict = [json.loads(shapely.to_geojson(geom))] # 提取指定范围内的数据 out_image, out_transform = mask(src, geom_dict, crop=True, all_touched=True) # 获取无数据值 nodata = src.nodata # 将三维数组转换为二维（如果是单波段） if out_image.ndim == 3: out_image = out_image[0] # 创建坐标网格 height, width = out_image.shape rows, cols = np.meshgrid(np.arange(height), np.arange(width), indexing='ij') # 计算每个像元的中心点坐标 lon_coords = [] lat_coords = [] for row in range(height): for col in range(width): # 计算像元中心点的坐标 lon, lat = rasterio.transform.xy(out_transform, row + 0.5, col + 0.5) lon_coords.append(lon) lat_coords.append(lat) # 准备数据 data = { '行号': rows.flatten(), '列号': cols.flatten(), '经度': lon_coords, '纬度': lat_coords, '栅格值': out_image.flatten() } # 创建DataFrame df = pd.DataFrame(data) # 过滤掉无数据值 if nodata is not None: df = df[df['栅格值'] != nodata] # 保存到Excel df.to_excel(output_excel_path, index=False, engine='openpyxl') print(f"成功提取数据到: {output_excel_path}") print(f"提取了 {len(df)} 个有效像元") print(f"数据范围: 经度 {df['经度'].min():.6f} - {df['经度'].max():.6f}") print(f" 纬度 {df['纬度'].min():.6f} - {df['纬度'].max():.6f}") return df except Exception as e: print(f"处理过程中出现错误: {e}") return None def batch_process_tif_files(input_folder, output_folder, min_lon, max_lon, min_lat, max_lat, resolution_meters=4000): """ 批量处理文件夹中的所有TIF文件 参数: input_folder: 输入文件夹路径 output_folder: 输出文件夹路径 min_lon, max_lon: 经度范围 min_lat, max_lat: 纬度范围 resolution_meters: 栅格分辨率 """ # 确保输出文件夹存在 if not os.path.exists(output_folder): os.makedirs(output_folder) print(f"创建输出文件夹: {output_folder}") # 查找所有的TIF文件 tif_files = glob.glob(os.path.join(input_folder, "*.tif")) print(f"找到 {len(tif_files)} 个TIF文件") # 处理每个文件 for tif_file in tif_files: # 获取文件名（不含扩展名） file_name = os.path.splitext(os.path.basename(tif_file))[0] # 构建输出Excel文件路径 output_excel = os.path.join(output_folder, f"{file_name}.xlsx") print(f"\n开始处理: {tif_file}") print(f"输出文件: {output_excel}") # 执行提取 result_df = extract_raster_to_excel( raster_path=tif_file, output_excel_path=output_excel, min_lon=min_lon, max_lon=max_lon, min_lat=min_lat, max_lat=max_lat, resolution_meters=resolution_meters ) if result_df is not None: print(f"成功处理: {file_name}") else: print(f"处理失败: {file_name}") # 使用示例 if __name__ == "__main__": # 设置输入和输出文件夹路径 input_folder = "D:/20240809example/CTT/" # 替换为您的TIF文件所在文件夹 output_folder = "D:/20240809example/CTT/提取结果/" # 输出Excel文件的文件夹 # 经纬度范围 min_longitude = 114.2692 max_longitude = 117.2414 min_latitude = 43.2313 max_latitude = 45.32772 # 执行批量处理 batch_process_tif_files( input_folder=input_folder, output_folder=output_folder, min_lon=min_longitude, max_lon=max_longitude, min_lat=min_latitude, max_lat=max_latitude, resolution_meters=4000 ) print("\n批量处理完成！")
浙公网安备 33010602011771号