建立在另一篇《python版GEE批量提取采样点最临近时间内landsat中心像元的反射率》的基础上来求的,由于上一篇已经求出来影像ID以及采样点中心像元最中心的坐标经纬度ImageCoords,根据像元之间的距离建立缓冲区圆,将圆变为矩形再变为数组,提取周围25个像元的反射率。
如图,5*5=25个像元(landsat一个像元分辨率为30m),知采样点中心像元最中心O点的坐标经纬度ImageCoords,建立一个缓冲区圆使其能够缓冲到附近的25个像元(不能超过25个像元要求建立的缓冲区圆半径小于OA距离,即小于75m;不能少于25个像元要求建立的缓冲区圆半径大于OB,即大于45倍根2;所以缓冲区半径要求大于等于64,小于75)
开始来写代码
1、导入必要的库
import ee
import geemap
import tqdm #显示进度条
import pandas as pd
import numpy as np
geemap.set_proxy(****)#端口号
ee.Initialize()
2*去云和转换计算,可要可不要,我写这段的时候没有要,因为上一篇笔记已经去云了,先记下来万一以后用得着
def maskL457sr(image): #定义一个影像去云,mask掉L457质量不好的函数
qaMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 2)).eq(0) #以QA波段为标准的mask,int('11111', 2)为31
saturationMask = image.select('QA_RADSAT').eq(0) #用于检查 QA_RADSAT 波段的每个像素值是否等于 0,0的部分为饱和,mask掉
# Apply the scaling factors to the appropriate bands.
opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2) #将常见的光谱波段进行转换计算
thermalBand = image.select('ST_B6').multiply(0.00341802).add(149.0) #热红外波段转换计算
# Replace the original bands with the scaled ones and apply the masks.
return image.addBands(opticalBands, None, True) \
.addBands(thermalBand, None, True) \
.updateMask(qaMask) \
.updateMask(saturationMask) #将换算过的SR-B。 热红外波段添加到影像中,更新mask后的QA_PIXEL,QA_RADSAT,返回影像
3、建立缓冲区,后面需要3*3,7*7也可以依次类推改一改参数就好了
df = pd.read_excel(r'C:\Users\带有影像ID和像元中心点经纬度的表格.xlsx')
bands = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7'] #将这几个波段放进字典bands中
datas = [] #准备一个表格存储数据
for index, row in tqdm.tqdm(df.iterrows()): #tqdm.tqdm()创建进度条对象,df.iterrows()逐行遍历表格,返回每一行索引index和当前行数据row
try:
# 获取image影像
image = ee.Image(row['ImageId'])
# 定义选择反射率波段
# 对反射率值进行缩放(官方提供)
image = image.select(bands).multiply(0.0000275).add(-0.2)
# 获取点坐标位置
point = ee.Geometry.Point([row['centrol_lon'], row['centrol_lat']])
# 以该点为中心生成5*5的矩阵区域
area_of_interest = point.buffer(64).bounds() #这个值不一定是64,64到74每一个数都行,一般来说,前提是这里提供的经纬度必须得是像元最中心
# 根据5*5矩形区域将该范围影像的值转为矩阵 图像对象转换为 NumPy 数组
im_data = geemap.ee_to_numpy(image,region=area_of_interest, scale=30)
# 验证一遍
if im_data.shape[0]*im_data.shape[1] == 25:
datas.append({"number":row['number'],"coordinates":(row['longitude'], row['latitude']), "matrix":im_data})
except Exception as e:
print(f"{row['number']} 点矩阵大小完全不符合要求")
# print(datas)
4、提取反射率到表格里去
# 创建一个空的 DataFrame 用于存储所有数据
all_data = []
for data in datas:
x, y = data['coordinates']
matrix = data['matrix']
number = data['number']
# print(f"Number: {number}, Coordinates: ({x}, {y}), Matrix: {matrix}")
# 添加当前点的数据
for i in range(5): # 遍历每一行
for j in range(5): # 遍历每一列
row = [number,x, y, i + 1, j + 1] # 添加点坐标和像素位置,从1开始
row.extend(matrix[i, j, :]) # 添加该像素的所有波段值,matrix为三位数组,第三围表示波段
all_data.append(row)
# 添加空行
all_data.append([''] * (5 + len(bands)))
# 创建 DataFrame
columns = ['number','Longitude', 'Latitude', 'Row', 'Col'] + bands
df = pd.DataFrame(all_data, columns=columns)
# 写入 Excel 文件
output_file = r'C:\Users\需要保存的表格.xlsx'
df.to_excel(output_file, index=False)
print(f"Data has been written to {output_file}")
读取原始数据的表格是这样的前面的经纬度是采样点经纬度,后面的经纬度是开头提到的上一篇笔记提的最中心的经纬度,也可以两篇笔记的代码结合一下
结果表格是这样的,有行列信息,能清楚指导25个像元位置