fix sentinel1 gcps

main.py

@@ -196,7 +196,7 @@ if __name__ == "__main__":
         # "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"
+        "output_path":r"D:\hunan\S1A_IW_GRDH_1SDV_20231014T102655_20231014T102720_050760_061DF1_1AE9.SAFE\test\test7.tif"
     }
 
     if input_config["satellite"] in visible_satellite:

utils/SarProcessor.py

@@ -4,6 +4,7 @@ from osgeo import gdal,osr
 from scipy.ndimage import uniform_filter, gaussian_filter
 from skimage.filters import threshold_otsu
 import os
+import xml.etree.ElementTree as ET
 # os.environ["PROJ_LIB"] = r"D:\Miniconda3\Library\share\proj\proj.db"
 
 class SarProcessor():
@@ -84,16 +85,41 @@ class SarProcessor():
         return water_mask
     
 
-    def mask_save(self,mask,output_file):
-        #防止ZIPDecode:Decoding 问题
-        self.profile.update(
-            {'compress':'LZW'}
-        )
-        with rasterio.open(output_file, 'w', **self.profile) as dst:
-            if mask.ndim != 2:
-                raise ValueError("Array must be 2D")
-            dst.write(mask, 1) 
-            print(self.profile)
+    # def mask_save(self,mask,output_file):
+    #     #防止ZIPDecode:Decoding 问题
+    #     self.profile.update(
+    #         {'compress':'LZW'}
+    #     )
+    #     with rasterio.open(output_file, 'w', **self.profile) as dst:
+    #         if mask.ndim != 2:
+    #             raise ValueError("Array must be 2D")
+    #         dst.write(mask, 1) 
+    #         print(self.profile)
+
+    def parse_kml_corners(self):
+        xml_path = os.path.join(os.path.dirname(os.path.dirname(self.vv_band_path)),'preview','map-overlay.kml')
+        # 加载并解析 KML 文件
+        print('xml_path',xml_path)
+        tree = ET.parse(xml_path)
+        root = tree.getroot()
+
+        # 查找 LatLonQuad 中的 coordinates 标签
+        coordinates_text = root.find('.//{*}LatLonQuad/{*}coordinates').text.strip()
+
+        # 分割出四个坐标点并转换为浮点数
+        coords = coordinates_text.split()
+        ul_lon, ul_lat = map(float, coords[0].split(','))  # 左上角
+        ur_lon, ur_lat = map(float, coords[1].split(','))  # 右上角
+        lr_lon, lr_lat = map(float, coords[2].split(','))  # 右下角
+        ll_lon, ll_lat = map(float, coords[3].split(','))  # 左下角
+
+        # 返回结果
+        return {
+            'ul': (ul_lon, ul_lat),
+            'ur': (ur_lon, ur_lat),
+            'lr': (lr_lon, lr_lat),
+            'll': (ll_lon, ll_lat)
+        }
 
 
     def write_tiff(self, mask, output_file):
@@ -101,10 +127,20 @@ class SarProcessor():
         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   # 左上角经度、纬度
+        coor_dict = self.parse_kml_corners()
+        ll_lon, ll_lat = coor_dict['ll'][0],coor_dict['ll'][1]   # 左下角经度、纬度
+        lr_lon, lr_lat = coor_dict['lr'][0],coor_dict['lr'][1]    # 右下角经度、纬度
+        ur_lon, ur_lat = coor_dict['ur'][0],coor_dict['ur'][1]    # 右上角经度、纬度
+        ul_lon, ul_lat = coor_dict['ul'][0],coor_dict['ul'][1]    # 左上角经度、纬度
+
+        print(ll_lon, ll_lat)
+        print(lr_lon, lr_lat)
+        print(ur_lon, ur_lat)
+        print(ul_lon, ul_lat)
+        # 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),  # 左下角