226 lines
9.6 KiB
Python
226 lines
9.6 KiB
Python
from osgeo import gdal
|
||
import cv2
|
||
import numpy as np
|
||
from osgeo import gdal
|
||
import ephem
|
||
|
||
# 水体掩膜
|
||
def is_water(reflectance_band4, reflectance_band5, threshold_reflectance_4=0.05, threshold_reflectance_5=0.15):
|
||
if reflectance_band4 < threshold_reflectance_4 and reflectance_band5 < threshold_reflectance_5:
|
||
return True
|
||
else:
|
||
return False
|
||
#云掩膜
|
||
def is_cloud(reflectance_band3, reflectance_band5, band14_00_00_temperature, threshold_reflectance_35=0.9, threshold_temperature_14=265):
|
||
if reflectance_band3 + reflectance_band5 > threshold_reflectance_35 and band14_00_00_temperature < threshold_temperature_14:
|
||
return True
|
||
else:
|
||
return False
|
||
# 太阳光耀
|
||
def calculate_solar_azimuth(latitude, longitude, observation_time):
|
||
observer = ephem.Observer()
|
||
observer.lat = str(latitude)
|
||
observer.lon = str(longitude)
|
||
observer.date = observation_time # Set the observation time
|
||
|
||
sun = ephem.Sun()
|
||
sun.compute(observer)
|
||
|
||
azimuth = float(sun.az) * 180 / ephem.pi
|
||
if azimuth < 0:
|
||
azimuth += 360
|
||
|
||
return azimuth
|
||
|
||
def main():
|
||
# 打开影像文件
|
||
image_00_00_path = "E:/Vscode_code/HS_H09_20230801_0000/HS_H09_20230801_0000.tif" # 替换为你的影像路径
|
||
dataset_00 = gdal.Open(image_00_00_path)
|
||
|
||
if dataset_00 is None:
|
||
print("HS_H09_20230801_0000.tif.")
|
||
return
|
||
|
||
# 水体掩膜
|
||
#获取波段4和波段5的表观反射率
|
||
band4_00_00_radiance = dataset_00.GetRasterBand(4).ReadAsArray()
|
||
band5_00_00_radiance = dataset_00.GetRasterBand(5).ReadAsArray()
|
||
|
||
# Specify your thresholds
|
||
threshold_reflectance_4 = 0.05
|
||
threshold_reflectance_5 = 0.15
|
||
|
||
#生成云掩膜
|
||
#获取波段3和波段5的表观反射率，波段14的亮温
|
||
band3_00_00_radiance = dataset_00.GetRasterBand(3).ReadAsArray()
|
||
band5_00_00_radiance = dataset_00.GetRasterBand(5).ReadAsArray()
|
||
band14_00_00_temperature = dataset_00.GetRasterBand(14).ReadAsArray()
|
||
|
||
# Specify your thresholds
|
||
threshold_reflectance_35 = 0.9
|
||
threshold_temperature_14 = 265
|
||
|
||
#计算太阳方位角
|
||
# Get the geotransformation information
|
||
geo_transform = dataset_00.GetGeoTransform()
|
||
|
||
# Get image size (number of rows and columns)
|
||
rows = dataset_00.RasterYSize
|
||
cols = dataset_00.RasterXSize
|
||
|
||
# Iterate through each pixel and convert row/col to geographic coordinates
|
||
for row in range(rows):
|
||
for col in range(cols):
|
||
x_geo = geo_transform[0] + col * geo_transform[1] + row * geo_transform[2]
|
||
y_geo = geo_transform[3] + col * geo_transform[4] + row * geo_transform[5]
|
||
|
||
# Replace with the actual latitude, longitude, and observation time
|
||
latitude = y_geo
|
||
longitude = x_geo
|
||
observation_time = "2023/08/01/00/00/00" # Replace with the actual observation time
|
||
|
||
solar_azimuth = calculate_solar_azimuth(latitude, longitude, observation_time)
|
||
|
||
# Create an empty mask to store water,cloud,sun pixels
|
||
water_cloud_sun_mask_image = []
|
||
|
||
# Iterate through each pixel and determine if it's water,cloud,sun or not
|
||
for y in range(len(band3_00_00_radiance)):
|
||
row_mask = []
|
||
for x in range(len(band3_00_00_radiance[y])):
|
||
|
||
if is_water(band4_00_00_radiance[y, x], band5_00_00_radiance[y, x], threshold_reflectance_4, threshold_reflectance_5):
|
||
row_mask.append(0)
|
||
|
||
elif is_cloud(band3_00_00_radiance[y, x], band5_00_00_radiance[y, x], band14_00_00_temperature[y, x], threshold_reflectance_35, threshold_temperature_14):
|
||
row_mask.append(0)
|
||
|
||
elif 165 < solar_azimuth < 200:
|
||
row_mask.append(0)
|
||
|
||
else:
|
||
row_mask.append(255)
|
||
|
||
water_cloud_sun_mask_image.append(row_mask)
|
||
|
||
# 遥感数据，假设你已经有了一个4通道的温度数据和对应的14通道数据
|
||
band7_00_00_temperature = dataset_00.GetRasterBand(7).ReadAsArray()
|
||
band14_00_00_temperature = dataset_00.GetRasterBand(14).ReadAsArray()
|
||
|
||
# 定义阈值
|
||
day_threshold = 360 # 白天温度阈值
|
||
night_threshold = 320 # 夜晚温度阈值
|
||
background_expansion = 7 # 背景区域扩展步长
|
||
|
||
# 计算晴空植被像元数量
|
||
# 计算NDVI
|
||
band4_00_00_radiance = dataset_00.GetRasterBand(4).ReadAsArray()
|
||
band3_00_00_radiance = dataset_00.GetRasterBand(3).ReadAsArray()
|
||
|
||
# 禁用特定的警告
|
||
np.seterr(divide='ignore', invalid='ignore')
|
||
|
||
# 计算NDVI，将无效值设为0
|
||
ndvi = np.where((band4_00_00_radiance + band3_00_00_radiance) != 0, (band4_00_00_radiance - band3_00_00_radiance) / (band4_00_00_radiance + band3_00_00_radiance), 0)
|
||
|
||
# 遥感数据，假设你已经有了一个植被指数的数据，比如Normalized Difference Vegetation Index (NDVI)
|
||
ndvi_data = ndvi # 替换为你的植被指数数据
|
||
|
||
# 定义阈值，这个阈值通常需要根据具体的数据和应用场景来设置
|
||
vegetation_threshold = 0.2 # 例如，设置为0.2表示NDVI值大于0.2的像元被认为是植被
|
||
|
||
# 定义植被列表
|
||
vegetation_band7_00_00 = []
|
||
vegetation_band14_00_00 = []
|
||
|
||
# 定义疑似高温火点列表，初始化为白色 (255)
|
||
suspicious_hotspots = np.full_like(band7_00_00_temperature, 255, dtype=np.uint8)
|
||
|
||
# 遍历每个像元
|
||
for row in range(len(band7_00_00_temperature)):
|
||
for col in range(len(band7_00_00_temperature[0])):
|
||
if water_cloud_sun_mask_image[row][col] == 255:
|
||
|
||
# 获取当前像元的温度数据
|
||
band7_00_00_temperature_current = band7_00_00_temperature[row][col]
|
||
band14_00_00_temperature_current = band14_00_00_temperature[row][col]
|
||
|
||
# 检查白天和夜晚的温度阈值条件
|
||
if band7_00_00_temperature_current > day_threshold or band7_00_00_temperature_current > night_threshold :
|
||
|
||
# 初始化邻域扩展步长
|
||
expansion = background_expansion
|
||
|
||
# 初始化晴空植被像元数量
|
||
clear_vegetation_count = 0
|
||
|
||
# 开始扩展邻域，直到满足晴空植被像元数量条件
|
||
while clear_vegetation_count < (expansion * expansion * 0.2):
|
||
|
||
# 获取邻域像元数据
|
||
neighborhood_band7_00_00 = band7_00_00_temperature[
|
||
row - expansion // 2 : row + expansion // 2 + 1,
|
||
col - expansion // 2 : col + expansion // 2 + 1,
|
||
]
|
||
|
||
# 计算左上角和右上角的行列索引
|
||
top_left_row = max(row - expansion // 2, 0)
|
||
top_left_col = max(col - expansion // 2, 0)
|
||
bottom_right_row = min(row + expansion // 2, 6001 - 1)
|
||
bottom_right_col = min(col + expansion // 2 + 1, 6001)
|
||
|
||
# 遍历每个像元
|
||
for row in range(top_left_row, +bottom_right_row + 1):
|
||
for col in range(top_left_col, bottom_right_col):
|
||
# 获取当前像元的植被指数值
|
||
ndvi_value = ndvi_data[row][col]
|
||
|
||
# 判断是否为晴空植被像元
|
||
if ndvi_value > vegetation_threshold:
|
||
# 如果是晴空植被像元，将其添加到晴空植被像元数量中
|
||
clear_vegetation_count += 1
|
||
|
||
# 计算背景晴空植被像元band7和band14亮温
|
||
vegetation_band7_00_00.append(band7_00_00_temperature[row][col])
|
||
vegetation_band14_00_00.append(band14_00_00_temperature[row][col])
|
||
|
||
# 如果还不满足条件，增加邻域扩展步长
|
||
if clear_vegetation_count < (expansion * expansion * 0.2):
|
||
expansion += 2 # 可以按需调整步长，例如+2
|
||
|
||
# 如果扩展步长超过最大值（例如，19x19），则退出循环
|
||
if expansion > 19:
|
||
break
|
||
|
||
# 计算温度差和平均温度差
|
||
delta_B7_14 = band7_00_00_temperature_current - band14_00_00_temperature_current
|
||
delta_B7_11_bg = np.array(vegetation_band7_00_00) - np.array(vegetation_band14_00_00)
|
||
|
||
# 检查疑似高温火点条件
|
||
if band7_00_00_temperature_current > np.mean(vegetation_band7_00_00) + 10 and delta_B7_14 > np.mean(delta_B7_11_bg) + 8:
|
||
# 将当前像元标记为疑似高温火点
|
||
suspicious_hotspots[row][col] = np.array( 0, dtype=np.uint8) # 黑色
|
||
|
||
suspicious_hotspots= np.array(suspicious_hotspots, dtype=np.uint8)
|
||
|
||
# 保存图像到本地
|
||
cv2.imwrite("E:/Vscode_code/threshold_adapt_fireimage/suspicious_hotspots.tif", suspicious_hotspots)
|
||
|
||
# Open target TIFF for update
|
||
target_dataset = gdal.Open("E:/Vscode_code/threshold_adapt_fireimage/suspicious_hotspots.tif", gdal.GA_Update)
|
||
dataset_00_geo_transform = dataset_00.GetGeoTransform()
|
||
|
||
# Apply the source geo-transform to the target TIFF
|
||
target_dataset.SetGeoTransform(dataset_00_geo_transform)
|
||
|
||
|
||
# Close the datasets to save changes
|
||
dataset_00 = None
|
||
target_dataset = None
|
||
|
||
cv2.waitKey(0)
|
||
cv2.destroyAllWindows()
|
||
|
||
if __name__ == "__main__":
|
||
main()
|