基于矢量点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)