依据矢量点提取多张多波段影像的像素值，并将其存入EXCEL表格

最新推荐文章于 2024-09-09 00:00:00 发布

qq_43444176

最新推荐文章于 2024-09-09 00:00:00 发布

阅读量246

点赞数 1

文章标签： python

本文链接：https://blog.csdn.net/qq_43444176/article/details/130069917

版权

基于矢量点shp文件提取多张影像的对应点的像素值（多变量），并从影像文件名中获取日期-时间，将其与各变量一并存入EXCEL表格。

import rasterio
import geopandas as gpd
import numpy as np
import pandas as pd
import os
import glob
from osgeo import gdal
from datetime import datetime
from pathlib import Path
from datetime import datetime, timedelta
import geopandas as gpd
import rasterio

# Load the vector points as a GeoDataFrame
points = gpd.read_file(r'F:\Solar_Radiance\BSRN_SSR\Station_Point\Alberts_Point\54421new.shp')

# create an empty list to store the results
results = []
img_folder = r'F:\Solar_Radiance\BSRN_SSR\Factor'
for img_path in glob.glob(os.path.join(img_folder, '*.tif')):
    # Load the multi-band TIFF image
    with rasterio.open(img_path) as src:
        data = src.read()
        transform = src.transform
        crs = src.crs
        num_bands = src.count

    # Reproject the vector points to the same CRS as the TIFF image, if necessary
    if points.crs != crs:
        points = points.to_crs(crs)

    # Extract the pixel values from the TIFF image based on the vector points
    values = []

    #这里flag的设置和循环里的if都是针对仅包含单个点的shp文件，多点文件需要再考虑
    flag = 0
    for idx, point in points.iterrows():
        col, row = rasterio.transform.rowcol(transform, point['geometry'].x, point['geometry'].y)
        if row < 0 or row >= data.shape[1] or col < 0 or col >= data.shape[2]:
            print(f"Skipping point: ({col}, {row}) is outside the raster extent.")
            flag = 1
            break
        pixel_values = data[:, row, col]
        values.append(pixel_values.flatten())

    if flag == 1 :
        continue

    date_str = (Path(img_path).stem)[1:8]
    time_str = (Path(img_path).stem)[8:12]
    formatted_time_str = time_str[:2] + ":" + time_str[2:] + ":00"
    time_obj = datetime.strptime(formatted_time_str, '%H:%M:%S')  # convert string to datetime object
    new_time_obj = time_obj + timedelta(hours=8)  # add 8 hours to the datetime object
    new_time_str = new_time_obj.strftime('%H:%M:%S')  # convert datetime object back to string
    print(new_time_str)

    # date = datetime.strptime(date_str, '%Y%m%d')
    # doy = (date - datetime(date.year, 1, 1)).days + 1
    row = [date_str, new_time_str] + np.concatenate(values).tolist()

    # add the row to the results list
    results.append(row)
    # Add the pixel values as new columns to the GeoDataFrame
    #for band_idx in range(num_bands):
        #points[f'band{band_idx+1}_value'] = np.array(values)[:, band_idx]

    # Save the GeoDataFrame to a new shapefile
    #points.to_file(r'F:\Solar_Radiance\BSRN_SSR\Factor_try\points_with_values.shp')

Factor_name_List = ['ETR','Tr','Tum','To','Tw','Ta','Taa','Rs','WSA','BSA','Cld_Par_Eff_Rad','Cld_Water_Path','CF']
df = pd.DataFrame(results, columns=['Date', 'Time'] + Factor_name_List)
df.to_excel(r'F:\Solar_Radiance\BSRN_SSR\Factor\54221.xlsx', index=False)

# Create a pandas DataFrame with the pixel values for every band
#df = pd.DataFrame(values, columns=[f'band{band_idx+1}' for band_idx in range(num_bands)])

# Write the DataFrame to an Excel file
#df.to_excel(r'F:\Solar_Radiance\BSRN_SSR\Factor_try\output.xlsx', index=False)