add sentinel1 gcp

main.py

@@ -145,30 +145,31 @@ def visible_go(input_config):
     return 'success'
 
 
-# def sar_go(input_config):
+def sar_go(input_config):
 
-#     band_path_dic = search_band(input_config,ref_config)
-#     satellite = input_config['satellite']
-#     print(band_path_dic)
+    band_path_dic = search_band(input_config,ref_config)
+    satellite = input_config['satellite']
+    print(band_path_dic)
 
-#     #load tif数据
-#     index = input_config["index"]
-#     sar_processor = SarProcessor(band_path_dict = band_path_dic,index = index )
-#     vv_array, vh_array = sar_processor.load_bands()
+    #load tif数据
+    index = input_config["index"]
+    sar_processor = SarProcessor(band_path_dict = band_path_dic,index = index )
+    vv_array, vh_array = sar_processor.load_bands()
 
-#     #lee滤波
-#     filtered_vv_array = sar_processor.lee_filter(vv_array) 
-#     filtered_vh_array = sar_processor.lee_filter(vh_array) 
+    #lee滤波
+    filtered_vv_array = sar_processor.lee_filter(vv_array) 
+    filtered_vh_array = sar_processor.lee_filter(vh_array) 
 
-#     #vv/vh比值
-#     water_ratio = sar_processor.compute_ratio(filtered_vv_array, filtered_vh_array)
+    #vv/vh比值
+    water_ratio = sar_processor.compute_ratio(filtered_vv_array, filtered_vh_array)
 
-#     #Ostu阈值法提取水体
-#     water_mask = sar_processor.extract_water(water_ratio)
-#     output_path = input_config["output_path"]
-#     os.makedirs(os.path.dirname(input_config["output_path"]), exist_ok=True)
-#     sar_processor.mask_save(water_mask,output_path)
-#     return 'sar success'
+    #Ostu阈值法提取水体
+    water_mask = sar_processor.extract_water(water_ratio)
+    output_path = input_config["output_path"]
+    os.makedirs(os.path.dirname(input_config["output_path"]), exist_ok=True)
+    # sar_processor.mask_save(water_mask,output_path)
+    sar_processor.write_tiff(water_mask,output_path)
+    return 'sar success'
 
 
 if __name__ == "__main__":
@@ -181,28 +182,28 @@ if __name__ == "__main__":
     #     "output_path": r"D:\hunan\50RKQ_l2a\test2"
     # }
 
-    input_config = {
-        "index": "NDVI",
-        "satellite": "landsat",
-        "input_path": r"D:\hunan\landsat_data\LC08_L2SP_122043_20231227_20240104_02_T1\LC08_L2SP_122043_20231227_20240104_02_T1_MTL.xml",
-        "output_path": r"D:\hunan\landsat_data\LC08_L2SP_122043_20231227_20240104_02_T1\test1"
-    }
-
     # input_config = {
-    #     "index": "NDWI",
-    #     "satellite": "sentinel1",
-    #     # "input_path":r"D:\hunan\s1_20240116_aws\S1A_IW_GRDH_1SDV_20231031T103557_20231031T103622_051008_062673_CCD9\manifest.safe",
-    #     # "output_path": r"D:\hunan\s1_20240116_aws\S1A_IW_GRDH_1SDV_20231031T103557_20231031T103622_051008_062673_CCD9\test\test_water3_mask.tif"
-    #     "input_path":r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\manifest.safe",
-    #     # "output_path":r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\test\test2.tif"
-    #     "output_path":r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\test"
+    #     "index": "NDVI",
+    #     "satellite": "landsat",
+    #     "input_path": r"D:\hunan\landsat_data\LC08_L2SP_122043_20231227_20240104_02_T1\LC08_L2SP_122043_20231227_20240104_02_T1_MTL.xml",
+    #     "output_path": r"D:\hunan\landsat_data\LC08_L2SP_122043_20231227_20240104_02_T1\test1"
     # }
 
+    input_config = {
+        "index": "NDWI",
+        "satellite": "sentinel1",
+        # "input_path":r"D:\hunan\s1_20240116_aws\S1A_IW_GRDH_1SDV_20231031T103557_20231031T103622_051008_062673_CCD9\manifest.safe",
+        # "output_path": r"D:\hunan\s1_20240116_aws\S1A_IW_GRDH_1SDV_20231031T103557_20231031T103622_051008_062673_CCD9\test\test_water3_mask.tif"
+        "input_path":r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\manifest.safe",
+        # "output_path":r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\test\test2.tif"
+        "output_path":r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\test\test6.tif"
+    }
+
     if input_config["satellite"] in visible_satellite:
         result = visible_go(input_config)
-    # elif input_config["satellite"] in sar_satellite:
-        # result = sar_go(input_config)
-    # print(result)
+    elif input_config["satellite"] in sar_satellite:
+        result = sar_go(input_config)
+    print(result)
 
 
     

utils/SarProcessor.py

 import rasterio
 import numpy as np
-from osgeo import gdal
+from osgeo import gdal,osr
 from scipy.ndimage import uniform_filter, gaussian_filter
 from skimage.filters import threshold_otsu
 import os
@@ -18,6 +18,8 @@ class SarProcessor():
     def load_bands(self):
         with rasterio.open(self.vv_band_path,'r') as src:
             crs = rasterio.crs.CRS.from_epsg(4326)
+            self.width = src.width
+            self.height = src.height
             print('crs',crs)
             print('crs',crs)
             src.profile.update(crs = crs)
@@ -93,17 +95,40 @@ class SarProcessor():
             dst.write(mask, 1) 
             print(self.profile)
 
-    # def write_tiff(self, mask, output_file):
-    #     driver = gdal.GetDriverByName('Gtiff')
-    #     out_dataset = driver.Create(output_file, self.width, self.height, 1, gdal.GDT_Byte)
-    #     out_dataset.SetGeoTransform(self.geo_transform)
+
+    def write_tiff(self, mask, output_file):
+        print('开始写出')
+        driver = gdal.GetDriverByName('Gtiff')
+        out_dataset = driver.Create(output_file, self.width, self.height, 1, gdal.GDT_Byte)
+        # 指定四个角的经纬度坐标（左下、右下、右上、左上）
+        ll_lon, ll_lat = 113.458282,25.395218  # 左下角经度、纬度
+        lr_lon, lr_lat = 115.971077,25.815006    # 右下角经度、纬度
+        ur_lon, ur_lat = 115.672813,27.319731    # 右上角经度、纬度
+        ul_lon, ul_lat = 113.126534,26.902760   # 左上角经度、纬度
+
+        gcp_list = [
+            gdal.GCP(ll_lon, ll_lat, 0, 0, 0),  # 左下角
+            gdal.GCP(lr_lon, lr_lat, 0, self.width, 0),  # 右下角
+            gdal.GCP(ur_lon, ur_lat, 0, self.width, self.height),  # 右上角
+            gdal.GCP(ul_lon, ul_lat, 0, 0, self.height)  # 左上角
+        ]
+        # 应用 GCP 点
+        print('应用 GCP 点')
+        sr = osr.SpatialReference()
+        sr.SetWellKnownGeogCS('WGS84')
+        out_dataset.SetGCPs(gcp_list, sr.ExportToWkt())
+        # 根据 GCP 计算仿射变换矩阵
+        print('开始计算仿射变换矩阵')
+        self.geo_transform = gdal.GCPsToGeoTransform(gcp_list)
+        out_dataset.SetGeoTransform(self.geo_transform)
         # out_dataset.SetProjection(self.projection)
-    #     out_band = out_dataset.GetRasterBand(1)
-    #     out_band.WriteArray(mask)
-    #     out_band.SetNoDataValue(0)
-    #     out_band.FlushCache()
-    #     out_dataset = None
-    #     print('sar success')
+        out_band = out_dataset.GetRasterBand(1)
+        out_band.WriteArray(mask)
+        out_band.SetNoDataValue(0)
+        out_band.FlushCache()
+        out_dataset = None
+        print('sar success')
+
 
 if __name__ == "__main__":
     path = r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\manifest.safe"