From 1b28df6d76321ae1645c81bf50540c194bd2d236 Mon Sep 17 00:00:00 2001 From: Wc098 <1832724385@qq.com> Date: Sun, 26 Nov 2023 23:00:11 +0800 Subject: [PATCH] =?UTF-8?q?=E6=A3=AE=E6=9E=97=E7=81=AB=E7=82=B9=E8=AF=86?= =?UTF-8?q?=E5=88=AB=E7=AE=97=E6=B3=95=E4=BB=A3=E7=A0=81?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .vscode/launch.json | 12 +++ threshold.py | 74 +++++++++++++++ threshold_adapt.py | 225 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 311 insertions(+) create mode 100644 .vscode/launch.json create mode 100644 threshold.py create mode 100644 threshold_adapt.py diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..6cf1814 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,12 @@ +{ + "configurations": [ + + { + "name": "Python: File", + "type": "python", + "request": "launch", + "program": "${file}", + "justMyCode": true + } + ] +} \ No newline at end of file diff --git a/threshold.py b/threshold.py new file mode 100644 index 0000000..781c772 --- /dev/null +++ b/threshold.py @@ -0,0 +1,74 @@ +from osgeo import gdal +import cv2 +import numpy as np +from osgeo import gdal + +def main(): + # 打开影像文件 + image_00_00_path = "E:/vscode_code/fire_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 + + image_00_10_path = "E:/vscode_code/fire_code/HS_H09_20230801_0010/HS_H09_20230801_0010.tif" # 替换为你的影像路径 + dataset_10 = gdal.Open(image_00_10_path) + + if dataset_10 is None: + print("HS_H09_20230801_0010.tif.") + return + + # 获取00时00分 band 7 和 band 14 的 temperature 数据 + band7_00_00_temperature = dataset_00.GetRasterBand(7).ReadAsArray() + band14_00_00_temperature = dataset_00.GetRasterBand(14).ReadAsArray() + + # 获取00时10分 band 7 和 band 14 的 temperature 数据 + band7_00_10_temperature = dataset_10.GetRasterBand(7).ReadAsArray() + band14_00_10_temperature = dataset_10.GetRasterBand(14).ReadAsArray() + + #计算band7，band14在00分到10分的差值 + band7_temperature_difference = band7_00_10_temperature - band7_00_00_temperature + band14_temperature_difference = band14_00_10_temperature - band14_00_00_temperature + + #计算band7与band14在00分和10分的差值 + band7_band14_00temperature_difference = band7_00_00_temperature - band14_00_00_temperature + band7_band14_10temperature_difference = band7_00_10_temperature - band14_00_10_temperature + + #计算10分的band7与00分的band14的差值 + band7_band14_0010temperature_difference = band7_00_10_temperature - band14_00_00_temperature + + #判断火点依据 + fire_point_image = [] + + for y in range(len(band7_00_00_temperature)): + row = [] + for x in range(len(band7_00_00_temperature[y])): + if band7_00_00_temperature[y][x] > 260 and band7_00_10_temperature[y][x] > 260 and band7_temperature_difference[y][x]>15 and band7_band14_10temperature_difference[y][x] >10 and band14_temperature_difference[y][x] >-1 and band7_band14_0010temperature_difference[y][x] >12: + row.append(0) # 白色 + else: + row.append(225) # 黑色 + fire_point_image.append(row) + + fire_point_image = np.array(fire_point_image, dtype=np.uint8) + + # 保存火点图像到本地 + cv2.imwrite("fire_point_image_2.tif", fire_point_image) + + # Open target TIFF for update + target_dataset = gdal.Open("fire_point_image_2.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 + dataset_10 = None + target_dataset = None + + cv2.waitKey(0) + cv2.destroyAllWindows() + +if __name__ == "__main__": + main() diff --git a/threshold_adapt.py b/threshold_adapt.py new file mode 100644 index 0000000..1ff2543 --- /dev/null +++ b/threshold_adapt.py @@ -0,0 +1,225 @@ +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()