修改盐碱度分块后取值越界问题

Ortho/Ortho.xml

@@ -1,7 +1,7 @@
 <?xml version='1.0' encoding='utf-8'?>
 <Root>
   <TaskID>CSAR_202107275419_0001-0</TaskID>
-  <WorkSpace>E:\Result_GF3\</WorkSpace>
+  <WorkSpace>D:\micro\WorkSpace\</WorkSpace>
   <AlgCompt>
     <DataTransModel>File</DataTransModel>
     <Artificial>ElementAlg</Artificial>
@@ -45,7 +45,7 @@
         <ParaType>File</ParaType>
         <DataType>tar.gz</DataType>
         <ParaSource>Cal</ParaSource>
-        <ParaValue>E:\GF3Data\soilMoisture\GF3B_SYC_QPSI_008316_E116.1_N43.3_20230622_L1A_AHV_L10000202892.tar.gz</ParaValue>
+        <ParaValue>D:\micro\microproduct_depdence\GF3-Deformation\download\cls\GF3_SAY_FSI_001614_E113.2_N34.5_20161129_L1A_HHHV_L10002015686.tar.gz</ParaValue>
         <EnModification>True</EnModification> 
         <EnMultipleChoice>False</EnMultipleChoice> 
         <Control>File</Control> 
@@ -58,9 +58,9 @@
         <ParaChsName>DEM数字高程影像</ParaChsName>
         <Description>30m分辨率DEM数字高程影像tif</Description>
         <ParaType>File</ParaType>
-        <DataType>zip</DataType>
+        <DataType>File</DataType>
         <ParaSource>Cal</ParaSource>
-        <ParaValue>E:\GF3Data\ortho\CSAR_ortho_ASTGTM2_N42E115_dem.zip;E:\GF3Data\ortho\CSAR_ortho_ASTGTM2_N42E116_dem.zip;E:\GF3Data\ortho\CSAR_ortho_ASTGTM2_N43E115_dem.zip;E:\GF3Data\ortho\CSAR_ortho_ASTGTM2_N43E116_dem.zip</ParaValue>
+        <ParaValue>D:\micro\microproduct_depdence\GF3-Deformation\dem</ParaValue>
         <EnModification>True</EnModification> 
         <EnMultipleChoice>True</EnMultipleChoice> 
         <Control>File</Control> 
@@ -92,7 +92,7 @@
         <ParaType>File</ParaType>
         <DataType>tar.gz</DataType>
 		<ParaSource>Cal</ParaSource>
-        <ParaValue>E:\Result_GF3\Ortho\Output\GF3B_SYC_QPSI_008316_E116.1_N43.3_20230622_L1A_AHV_L10000202892-ortho.tar.gz</ParaValue>
+        <ParaValue>D:\micro\WorkSpace\Ortho\Output\GF3_SAY_FSI_001614_E113.2_N34.5_20161129_L1A_HHHV_L10002015686-ortho.tar.gz</ParaValue>
         <MaxValue>DEFAULT</MaxValue>
         <MinValue>DEFAULT</MinValue>
         <OptionValue>DEFAULT</OptionValue>

dem-C-SAR/Dem.xml

@@ -71,7 +71,7 @@
                 <ParaType>value</ParaType>
                 <DataType>string</DataType>
                 <ParaSource>Man</ParaSource>
-                <ParaValue>34.60;34.67;113.05;113.18</ParaValue>
+                <ParaValue>34.6;34.63;113.1;113.15</ParaValue>
                 <EnModification>True</EnModification>
                 <EnMultipleChoice>True</EnMultipleChoice>
                 <Control>UploadInput</Control>

dem-C-SAR/DemMain.py

@@ -458,6 +458,7 @@ class DemMain:
 
         # os.chdir(isce_exe_dir)
         cmd = "stackStripMap.exe -s {} -w {} -d {} -m {} -a {} -r {} -x {} -u 'snaphu' --nofocus".format(out_slc_dir, isce_work_space, dem_path, main_img, 3, 3, box)
+        # cmd = "stackStripMap.exe -s {} -w {} -d {} -m {} -a {} -r {} -u 'snaphu' --nofocus".format(out_slc_dir, isce_work_space, dem_path, main_img, 3, 3)
         logger.info('stackStripMap_cmd:{}'.format(cmd))
         result = os.system(cmd)
         logger.info('cmd_result:{}'.format(result))

dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Orbit/Orbit.py

@@ -234,6 +234,15 @@ class Orbit(Component):
         self._orbitSource = source or None
         self._referenceFrame = None
         self._stateVectors.configure()
+        """ 自定义 """
+        self.oribitStartTime = None
+        self.A_arr = None
+        self.polynum = None
+        """ 多普勒参数"""
+        self.refrenceTime = None
+        self.dopperPoly = []
+        """ 卫星名称"""
+        self.sessionMode = None
         # self._stateVectors = stateVectors or []
         self._cpStateVectors = []
         type(self._stateVectors)
@@ -260,6 +269,69 @@ class Orbit(Component):
     #    #pack the orbit into stateVectors
     #    self._packOrbit(cpStateVectors[0], cpStateVectors[1], cpStateVectors[2], cpStateVectors[3])
 
+    def setsessionMode(self, SessionMode):
+        self.sessionMode = SessionMode
+
+    def setPolyParams(self, polynum, orbitStartTime, A_arr):
+        self.polynum = polynum
+        self.oribitStartTime = orbitStartTime
+        self.A_arr = A_arr
+
+    def setDoppler(self, refrenceTime, dopperPoly):
+        self.refrenceTime = refrenceTime
+        self.dopperPoly = dopperPoly
+        pass
+
+    def getSatelliteSpaceState(self, time_float_datetime):
+        '''
+        逐像素求解
+        根据时间戳，返回对应时间的卫星的轨迹状态，会自动计算与起算时间之差
+        args:
+            time_float:时间戳
+        return：
+            State_list:[time,Xp,Yp,Zp,Vx,Vy,Vz]
+        '''
+        time_float = time_float_datetime.timestamp()
+        time_float = time_float - self.oribitStartTime
+        #
+        px = 0
+        py = 0
+        pz = 0
+        vx = 0
+        vy = 0
+        vz = 0
+        for ii in range(self.polynum):
+            px += self.A_arr[ii][0] * time_float ** ii
+            py += self.A_arr[ii][1] * time_float ** ii
+            pz += self.A_arr[ii][2] * time_float ** ii
+            vx += self.A_arr[ii][3] * time_float ** ii
+            vy += self.A_arr[ii][4] * time_float ** ii
+            vz += self.A_arr[ii][5] * time_float ** ii
+        sv = StateVector()
+        sv.setTime(time_float_datetime)
+        sv.setPosition([px, py, pz])
+        sv.setVelocity([vx, vy, vz])
+        return sv
+
+    def getTimeOrbits(self, UTCStartTime, UTCEndTime, orbitnum=1000):
+        #
+        startTime_stamp = datetime.datetime.strptime(UTCStartTime, "%Y-%m-%dT%H:%M:%S.%f") - datetime.timedelta(
+            seconds=30)
+        endTime_stamp = datetime.datetime.strptime(UTCEndTime, "%Y-%m-%dT%H:%M:%S.%f") + datetime.timedelta(seconds=30)
+        if startTime_stamp > endTime_stamp:
+            raise
+        delta_t = (endTime_stamp.timestamp() - startTime_stamp.timestamp()) / orbitnum
+        extractOrbits = []
+        #
+        temptime = startTime_stamp
+        while temptime < endTime_stamp:
+            temptime = temptime + datetime.timedelta(seconds=delta_t)
+            newOrbit = self.getSatelliteSpaceState(temptime)
+            extractOrbits.append(newOrbit)
+        newOrbit = self.getSatelliteSpaceState(endTime_stamp)
+        extractOrbits.append(newOrbit)
+        return extractOrbits  # 扩展的轨道节点
+
     def adaptToRender(self):
         import copy
         # make a copy of the stateVectors to restore it after dumping
@@ -376,10 +448,11 @@ class Orbit(Component):
         @raises NotImplementedError: if the desired interpolation method
         cannot be decoded
         """
+        if self.sessionMode is None:
             if time not in self:
                 raise ValueError(
-                "Time stamp (%s) falls outside of the interpolation interval [%s:%s]" %
+                    "Time stamp (%s) falls outside of the interpolation interval [%s:%s]" % (
-                (time,self.minTime,self.maxTime)
+                    time, self.minTime, self.maxTime)
                 )
 
             if method == 'linear':
@@ -393,9 +466,13 @@ class Orbit(Component):
                     "Orbit interpolation type %s, is not implemented" % method
                 )
             return newSV
+        elif self.sessionMode == "GF3" or self.sessionMode == "GF3B" or self.sessionMode == "GF3C":
+            print("GF3轨道插值")
+            return self.getSatelliteSpaceState(time)
 
     ## Isn't orbit redundant? -compute the method based on name
     def interpolate(self, time, method='linear'):
+        if self.sessionMode is None:
             if time not in self:
                 raise ValueError("Time stamp (%s) falls outside of the interpolation interval [%s:%s]"
                                  % (time, self.minTime, self.maxTime))
@@ -406,8 +483,11 @@ class Orbit(Component):
             raise NotImplementedError(
                 "Orbit interpolation type %s, is not implemented" % method
             )
+        elif self.sessionMode == "GF3" or self.sessionMode == "GF3B" or self.sessionMode == "GF3C":
+            print("GF3轨道插值")
+            return self.getSatelliteSpaceState(time)
 
-    interpolateOrbit = interpolate
+    # interpolateOrbit = interpolate
 
     def _linearOrbitInterpolation(self,time):
         """
@@ -998,7 +1078,7 @@ class Orbit(Component):
     pointOnGround = rdr2geo
 
     def geo2rdr(self, llh, side=-1, planet=None,
-            doppler=None, wvl=None):
+                doppler=None, wvl=None, isLT1AB=True):
         '''
         Takes a lat, lon, height triplet and returns azimuth time and range.
         Assumes zero doppler for now.
@@ -1014,48 +1094,180 @@ class Orbit(Component):
             refElp = Planet(pname='Earth').ellipsoid
         else:
             refElp = planet.ellipsoid
-
+        # print('llh', llh)
         xyz = refElp.llh_to_xyz(llh)
-
         delta = (self.maxTime - self.minTime).total_seconds() * 0.5
-        tguess = self.minTime + datetime.timedelta(seconds = delta)
+        tguess = self.minTime  # + datetime.timedelta(seconds = delta)
+        # print("Orbits.py 1024-----------------------------------------------")
+        # print("self.maxTime", self.maxTime)
+        # print("self.minTime", self.minTime)
+        # print(delta)
+        # print(tguess)
+
+        LIGHTSPEED = 299792458
+        if wvl == 0:
+            isLT1AB = False
+        if isLT1AB and (self.sessionMode == "GF3" or self.sessionMode == "GF3B" or self.sessionMode == "GF3C"):  # 专门针对 LT1AB
+            print("LT1AB orbit.....")
+            dt = 0.0001
+            outOfBounds = False
+            for ii in range(51):
+                try:
+                    sv = self.getSatelliteSpaceState(tguess + datetime.timedelta(seconds=dt))  # 获取卫星的 位置、速度
+                except Exception as e:
+                    print(e)
+                    outOfBounds = True
+                    break
+
+                pos1 = np.array(sv.getPosition())  # 卫星坐标
+                vel1 = np.array(sv.getVelocity())  # 卫星速度
+                dr1 = pos1 - xyz
+                rng1 = np.linalg.norm(dr1)  # 斜距
+
+                # ((R_s1.array() * V_s1.array()).rowwise().sum().array() * (-2) / (R * this->lamda))[0];
+                FdTheory1 = -2 / (rng1 * wvl) * np.dot(dr1, vel1)
+
+                try:
+                    sv = self.getSatelliteSpaceState(tguess)
+                except Exception as e:
+                    print(e)
+                    outOfBounds = True
+                    break
+                pos2 = np.array(sv.getPosition())  # 卫星坐标
+                vel2 = np.array(sv.getVelocity())  # 卫星速度
+                dr2 = pos2 - xyz
+                rng = np.linalg.norm(dr2)  # 斜距
+                FdTheory2 = -2 / (rng * wvl) * np.dot(dr2, vel2)
+                TSR = rng * 2 / LIGHTSPEED - self.refrenceTime  # nx1
+
+                FdNumerical = 0
+                # FdNumerical=FdNumerical+self.dopperPoly[0]*TSR**0
+                # FdNumerical=FdNumerical+self.dopperPoly[1]*TSR**1
+                # FdNumerical=FdNumerical+self.dopperPoly[2]*TSR**2
+                # FdNumerical=FdNumerical+self.dopperPoly[3]*TSR**3
+
+                fdopper_grad = (FdTheory1 - FdTheory2) / dt
+                inc_t = (FdTheory2 - FdNumerical) / fdopper_grad
+                # print(inc_t,rng,FdNumerical,FdTheory2,tguess,pos2)
+                tguess = tguess - datetime.timedelta(seconds=inc_t)
+
+                if abs(inc_t) < 1e-6:
+                    break
+                else:
+                    t_prev_guess = tguess
+            # print(outOfBounds)
+            # print("end ------------------------------------------------------------\n")
+            if outOfBounds:
+                raise Exception('Interpolation time out of bounds')
+        else:
             outOfBounds = False
-        # Start the previous guess tracking with dummy value
-        t_prev_guess = tguess + datetime.timedelta(seconds=10)
             for ii in range(51):
                 try:
                     sv = self.interpolateOrbit(tguess, method='hermite')
-            except:
+                except Exception as e:
+                    if self.sessionMode == "LT1A" or self.sessionMode == "LT1B":
+                        sv = self.getSatelliteSpaceState(tguess)  # 获取卫星的 位置、速度
+                    print(e)
                     outOfBounds = True
                     break
 
                 pos = np.array(sv.getPosition())
                 vel = np.array(sv.getVelocity())
 
+                # print("xyz", xyz)
+                # print("pos", pos)
                 dr = xyz - pos
-            rng = np.linalg.norm(dr)
+                rng = np.linalg.norm(dr)  # 求斜距
+                # print("rng", rng)
+
+                dopfact = np.dot(dr, vel)  # fd 公式
+                # print("dopfact", dopfact)
 
-            dopfact = np.dot(dr,vel)
                 fdop = doppler(DTU.seconds_since_midnight(tguess), rng) * wvl * 0.5
-            fdopder = (0.5*wvl*doppler(DTU.seconds_since_midnight(tguess),rng+10.0) - fdop) / 10.0
+                # print("doppler", doppler(DTU.seconds_since_midnight(tguess),rng))
+                # print("wvl", wvl)
+                print("fdop", fdop)
 
+                fdopder = (0.5 * wvl * doppler(DTU.seconds_since_midnight(tguess), rng + 10.0) - fdop) / 10.0
+                # print("doppler2", doppler(DTU.seconds_since_midnight(tguess),rng+10.0))
+                print("fdopder", fdopder)
                 fn = dopfact - fdop * rng
                 c1 = -np.dot(vel, vel)
                 c2 = (fdop / rng + fdopder)
-
+                # print("c1", c1)
+                # print("c2", c2)
                 fnprime = c1 + c2 * dopfact
-
+                # print("时间为", fn/fnprime)
+                # if abs(fn/fnprime) > 1:
+                #     break
                 tguess = tguess - datetime.timedelta(seconds=fn / fnprime)
-            if abs(tguess - t_prev_guess).total_seconds() < 5e-9:
+                # print("输出的tguess", tguess)
-                break
+            # print(outOfBounds)
-            else:
+            print("end ------------------------------------------------------------\n")
-                t_prev_guess = tguess
-
             if outOfBounds:
                 raise Exception('Interpolation time out of bounds')
 
         return tguess, rng
 
+    # def geo2rdr(self, llh, side=-1, planet=None,
+    #         doppler=None, wvl=None):
+    #     '''
+    #     Takes a lat, lon, height triplet and returns azimuth time and range.
+    #     Assumes zero doppler for now.
+    #     '''
+    #     from isceobj.Planet.Planet import Planet
+    #     from isceobj.Util.Poly2D import Poly2D
+    #     if doppler is None:
+    #         doppler = Poly2D()
+    #         doppler.initPoly(azimuthOrder=0, rangeOrder=0, coeffs=[[0.]])
+    #         wvl = 0.0
+    #
+    #     if planet is None:
+    #         refElp = Planet(pname='Earth'). ellipsoid
+    #     else:
+    #         refElp = planet.ellipsoid
+    #
+    #     xyz = refElp.llh_to_xyz(llh)
+    #
+    #     delta = (self.maxTime - self.minTime).total_seconds() * 0.5
+    #     tguess = self.minTime + datetime.timedelta(seconds = delta)
+    #     outOfBounds = False
+    #     # Start the previous guess tracking with dummy value
+    #     t_prev_guess = tguess + datetime.timedelta(seconds=10)
+    #     for ii in range(51):
+    #         try:
+    #             sv = self.interpolateOrbit(tguess, method='hermite')
+    #         except:
+    #             outOfBounds = True
+    #             break
+    #
+    #         pos = np.array(sv.getPosition())
+    #         vel = np.array(sv.getVelocity())
+    #
+    #         dr = xyz-pos
+    #         rng = np.linalg.norm(dr)
+    #
+    #         dopfact = np.dot(dr,vel)
+    #         fdop = doppler(DTU.seconds_since_midnight(tguess),rng) * wvl * 0.5
+    #         fdopder = (0.5*wvl*doppler(DTU.seconds_since_midnight(tguess),rng+10.0) - fdop) / 10.0
+    #
+    #         fn = dopfact - fdop * rng
+    #         c1 = -np.dot(vel, vel)
+    #         c2 = (fdop/rng + fdopder)
+    #
+    #         fnprime = c1 + c2 * dopfact
+    #
+    #         tguess = tguess - datetime.timedelta(seconds = fn/fnprime)
+    #         if abs(tguess - t_prev_guess).total_seconds() < 5e-9:
+    #             break
+    #         else:
+    #             t_prev_guess = tguess
+    #
+    #     if outOfBounds:
+    #         raise Exception('Interpolation time out of bounds')
+    #
+    #     return tguess, rng
+
 
     def exportToC(self, reference=None):
         from isceobj.Util import combinedlibmodule

dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Sensor/GF3_SLC.py

@@ -164,6 +164,160 @@ class SatelliteOrbit(object):
         return None
 
 
+# class SatelliteOrbitFitPoly(SatelliteOrbit):
+#     '''
+#     继承于SatelliteOribit类，为拟合多项式实现方法
+#     '''
+#
+#     def __init__(self) -> None:
+#         super().__init__()
+#         self.modelName="多项式"
+#         self.polynum=4
+#
+#     def ReconstructionSatelliteOrbit(self, GPSPoints_list, starttime):
+#         if len(GPSPoints_list)==2:
+#             self.polynum=1
+#             self.starttime = starttime
+#
+#             record_count = len(GPSPoints_list)
+#             time_arr = np.zeros((record_count, 1), dtype=np.float64)  # 使用np.float64只是为了精度高些；如果32位也能满足需求，请用32位
+#             state_arr = np.zeros((record_count, 6), dtype=np.float64)
+#             A_arr = np.zeros((self.polynum+1, 6), dtype=np.float64) # 四次项
+#             X=np.ones((record_count,self.polynum+1),dtype=np.float64) # 记录时间坐标
+#             # 将点记录转换为自变量矩阵、因变量矩阵
+#
+#             for i in range(record_count):
+#                 GPSPoint = GPSPoints_list[i]
+#                 time_ = GPSPoint[0] - self.starttime # 为了保证精度，对时间进行缩放
+#                 X[i,:]=np.array([1,time_])
+#                 state_arr[i, :] = np.array(GPSPoint[1:],dtype=np.float64).reshape(1,6) # 空间坐标
+#             self.model_f=[]
+#             for i in range(6):
+#                 Y = state_arr[:, i].reshape(-1,1)
+#                 A_arr[:,i]=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)[:,0]
+#
+#             self.A_arr=copy.deepcopy(A_arr.copy())
+#             return self.A_arr
+#         elif len(GPSPoints_list) > 6:
+#             self.polynum=4
+#             # 多项式的节点数，理论上是超过5个可以起算，这里为了精度选择10个点起算。
+#             # 多项式 XA=Y ==> A=(X'X)^X'Y，其中 A 为待求系数，X为变量，Y为因变量
+#             # 这里使用三次项多项式，共有6组参数。
+#             # 声明自变量，因变量，系数矩阵
+#             self.starttime = starttime
+#
+#             record_count = len(GPSPoints_list)
+#             time_arr = np.zeros((record_count, 1), dtype=np.float64)  # 使用np.float64只是为了精度高些；如果32位也能满足需求，请用32位
+#             state_arr = np.zeros((record_count, 6), dtype=np.float64)
+#             A_arr = np.zeros((self.polynum+1, 6), dtype=np.float64) # 四次项
+#             X=np.ones((record_count,self.polynum+1),dtype=np.float64) # 记录时间坐标
+#             # 将点记录转换为自变量矩阵、因变量矩阵
+#
+#             for i in range(record_count):
+#                 GPSPoint = GPSPoints_list[i]
+#                 time_ = GPSPoint[0] - self.starttime # 为了保证精度，对时间进行缩放
+#                 X[i,:]=np.array([1,time_,time_**2,time_**3,time_**4])
+#                 state_arr[i, :] = np.array(GPSPoint[1:],dtype=np.float64).reshape(1,6) # 空间坐标
+#             self.model_f=[]
+#             for i in range(6):
+#                 Y = state_arr[:, i].reshape(-1,1)
+#                 A_arr[:,i]=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)[:,0]
+#
+#             self.A_arr=copy.deepcopy(A_arr.copy())
+#             ''' 测试误差
+#             from matplotlib import pyplot
+#             label_list=['x','y','z','vx','vy','vz']
+#             color_list=['r','g','b','gold','gray','pink']
+#             pyplot.figure()
+#             for i in range(6):
+#                 Y = state_arr[:, i]
+#                 Y_predict=self.model_f[i](X)
+#                 pyplot.subplot(int("23{}".format(i+1)))
+#                 d=Y-Y_predict
+#                 pyplot.plot(X,d,label=label_list[i],color=color_list[i])
+#                 pyplot.title("max:{}".format(np.max(d)))
+#                 #self.model_f.append(interpolate.interp1d(X,Y,kind='cubic',fill_value='extrapolate'))
+#             pyplot.legend()
+#             pyplot.show()
+#             '''
+#             return self.A_arr
+#         else:
+#             self.A_arr = None
+#             return None
+#
+#     def SatelliteSpaceState(self, time_float):
+#         '''
+#         逐像素求解
+#         根据时间戳，返回对应时间的卫星的轨迹状态，会自动计算与起算时间之差
+#         args:
+#             time_float:时间戳
+#         return：
+#             State_list:[time,Xp,Yp,Zp,Vx,Vy,Vz]
+#         '''
+#         if self.model_f is None:
+#             return None
+#
+#         result_arr=np.zeros((1,7))
+#
+#         time_float = time_float - self.starttime
+#         result_arr[0,0]=time_float
+#         #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4
+#         time_float=np.array([1,time_float,time_float**2,time_float**3,time_float**4]).reshape(1,5)
+#         result_arr=np.matmul(time_float,self.A_arr)
+#         return [time_float,result_arr]
+#
+#     def getSatelliteSpaceState(self, time_array):
+#         '''
+#         矩阵求解
+#         根据时间戳矩阵，返回对应时刻的卫星空间状态（位置，速度），且会自动计算与起算时间之差
+#         args:
+#             time_array:nparray nx1 时间戳
+#         return:
+#             SatellitSpaceStateArray:nparray nx6 状态信息
+#         '''
+#         if self.model_f is None:
+#             return None  # 返回None,表示没有结果
+#         if self.polynum==4:
+#             n=time_array.shape[0]
+#             result_arr_=np.zeros((n,6),dtype=np.float64)
+#             time_float = time_array - self.starttime
+#             time_float=time_float.reshape(-1) # nx1
+#             time_arr=np.ones((time_float.shape[0],5)) # nx5
+#             time_arr[:,1]=time_float
+#             time_arr[:,2]=time_float**2
+#             time_arr[:,3]=time_float**3
+#             time_arr[:,4]=time_float**4
+#             result_arr_=np.matmul(time_arr,self.A_arr) # nx5  5x6
+#             #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4
+#             #result_arr=result_arr_
+#             return result_arr_ # 位置矩阵
+#         else:
+#             n=time_array.shape[0]
+#             result_arr_=np.zeros((n,6),dtype=np.float64)
+#             time_float = time_array - self.starttime
+#             time_float=time_float.reshape(-1) # nx1
+#             time_arr=np.ones((time_float.shape[0],self.polynum+1)) # nx5
+#             time_arr[:,1]=time_float
+#             result_arr_=np.matmul(time_arr,self.A_arr) # nx5  5x6
+#             #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4
+#             #result_arr=result_arr_
+#             return result_arr_ # 位置矩阵
+
+
+# def ReconstructionSatelliteOrbit(GPSPoints_list, starttime):
+#     '''
+#     构建卫星轨道
+#     args:
+#         GPSPoints_list:卫星轨道点
+#         starttime:起算时间
+#     '''
+#
+#     SatelliteOrbitModel = SatelliteOrbitFitPoly()
+#     if SatelliteOrbitModel.ReconstructionSatelliteOrbit(GPSPoints_list, starttime=starttime) is None:
+#         return None
+#     return SatelliteOrbitModel
+
+
 class SatelliteOrbitFitPoly(SatelliteOrbit):
     '''
     继承于SatelliteOribit类，为拟合多项式实现方法
@@ -199,7 +353,6 @@ class SatelliteOrbitFitPoly(SatelliteOrbit):
             self.A_arr = copy.deepcopy(A_arr.copy())
             return self.A_arr
         elif len(GPSPoints_list) > 6:
-            self.polynum=4
             # 多项式的节点数，理论上是超过5个可以起算，这里为了精度选择10个点起算。
             # 多项式 XA=Y ==> A=(X'X)^X'Y，其中 A 为待求系数，X为变量，Y为因变量
             # 这里使用三次项多项式，共有6组参数。
@@ -209,37 +362,21 @@ class SatelliteOrbitFitPoly(SatelliteOrbit):
             record_count = len(GPSPoints_list)
             time_arr = np.zeros((record_count, 1), dtype=np.float64)  # 使用np.float64只是为了精度高些；如果32位也能满足需求，请用32位
             state_arr = np.zeros((record_count, 6), dtype=np.float64)
-            A_arr = np.zeros((self.polynum+1, 6), dtype=np.float64) # 四次项
+            A_arr = np.zeros((self.polynum, 6), dtype=np.float64)  # 四次项
-            X=np.ones((record_count,self.polynum+1),dtype=np.float64) # 记录时间坐标
+            X = np.ones((record_count, self.polynum), dtype=np.float64)  # 记录时间坐标
             # 将点记录转换为自变量矩阵、因变量矩阵
 
             for i in range(record_count):
                 GPSPoint = GPSPoints_list[i]
                 time_ = GPSPoint[0] - self.starttime  # 为了保证精度，对时间进行缩放
-                X[i,:]=np.array([1,time_,time_**2,time_**3,time_**4])
+                X[i, :] = np.array(list(map(lambda ii: time_ ** ii, range(self.polynum))))
                 state_arr[i, :] = np.array(GPSPoint[1:], dtype=np.float64).reshape(1, 6)  # 空间坐标
             self.model_f = []
             for i in range(6):
                 Y = state_arr[:, i].reshape(-1, 1)
                 A_arr[:, i] = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T, X)), X.T), Y)[:, 0]
 
-            self.A_arr=copy.deepcopy(A_arr.copy())
+            self.A_arr = A_arr.copy()
-            ''' 测试误差            
-            from matplotlib import pyplot
-            label_list=['x','y','z','vx','vy','vz']
-            color_list=['r','g','b','gold','gray','pink']
-            pyplot.figure()
-            for i in range(6):
-                Y = state_arr[:, i]
-                Y_predict=self.model_f[i](X)
-                pyplot.subplot(int("23{}".format(i+1)))
-                d=Y-Y_predict
-                pyplot.plot(X,d,label=label_list[i],color=color_list[i])
-                pyplot.title("max:{}".format(np.max(d)))
-                #self.model_f.append(interpolate.interp1d(X,Y,kind='cubic',fill_value='extrapolate'))
-            pyplot.legend()
-            pyplot.show()
-            '''
             return self.A_arr
         else:
             self.A_arr = None
@@ -260,63 +397,88 @@ class SatelliteOrbitFitPoly(SatelliteOrbit):
         result_arr = np.zeros((1, 7))
 
         time_float = time_float - self.starttime
-        result_arr[0,0]=time_float
-        #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4
-        time_float=np.array([1,time_float,time_float**2,time_float**3,time_float**4]).reshape(1,5)
-        result_arr=np.matmul(time_float,self.A_arr)
-        return [time_float,result_arr]
 
-    def getSatelliteSpaceState(self, time_array):
+        #
-        '''
+        px = 0
-        矩阵求解
+        py = 0
-        根据时间戳矩阵，返回对应时刻的卫星空间状态（位置，速度），且会自动计算与起算时间之差
+        pz = 0
-        args:
+        vx = 0
-            time_array:nparray nx1 时间戳
+        vy = 0
-        return:
+        vz = 0
-            SatellitSpaceStateArray:nparray nx6 状态信息
+        for ii in range(self.polynum):
-        '''
+            px += self.A_arr[ii, 0] * time_float ** ii
-        if self.model_f is None:
+            py += self.A_arr[ii, 1] * time_float ** ii
-            return None  # 返回None,表示没有结果
+            pz += self.A_arr[ii, 2] * time_float ** ii
-        if self.polynum==4:
+            vx += self.A_arr[ii, 3] * time_float ** ii
-            n=time_array.shape[0]
+            vy += self.A_arr[ii, 4] * time_float ** ii
-            result_arr_=np.zeros((n,6),dtype=np.float64)
+            vz += self.A_arr[ii, 5] * time_float ** ii
-            time_float = time_array - self.starttime
+
-            time_float=time_float.reshape(-1) # nx1
+        return [time_float, [px, py, pz, vx, vy, vz]]
-            time_arr=np.ones((time_float.shape[0],5)) # nx5
+
-            time_arr[:,1]=time_float
+    def getTimeOrbitStamp(self, UTCStartTime_float):
-            time_arr[:,2]=time_float**2
+        sv = _StateVector()
-            time_arr[:,3]=time_float**3
+        temp_sv = self.SatelliteSpaceState(UTCStartTime_float)
-            time_arr[:,4]=time_float**4
+        sv.timeStamp = datetime.datetime.fromtimestamp(UTCStartTime_float)
-            result_arr_=np.matmul(time_arr,self.A_arr) # nx5  5x6
+        sv.xPosition = temp_sv[1][0, 0]
-            #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4
+        sv.yPosition = temp_sv[1][0, 1]
-            #result_arr=result_arr_
+        sv.zPosition = temp_sv[1][0, 2]
-            return result_arr_ # 位置矩阵
+        sv.xVelocity = temp_sv[1][0, 3]
-        else:
+        sv.yVelocity = temp_sv[1][0, 4]
-            n=time_array.shape[0]
+        sv.zVelocity = temp_sv[1][0, 5]
-            result_arr_=np.zeros((n,6),dtype=np.float64)
+        return sv
-            time_float = time_array - self.starttime
+
-            time_float=time_float.reshape(-1) # nx1
+    def getTimeOrbits(self, UTCStartTime, UTCEndTime, orbitnum=1000):
-            time_arr=np.ones((time_float.shape[0],self.polynum+1)) # nx5
+        #
-            time_arr[:,1]=time_float
+        startTime_stamp = datetime.datetime.strptime(UTCStartTime, "%Y-%m-%dT%H:%M:%S.%f").timestamp() - 1
-            result_arr_=np.matmul(time_arr,self.A_arr) # nx5  5x6
+        endTime_stamp = datetime.datetime.strptime(UTCEndTime, "%Y-%m-%dT%H:%M:%S.%f").timestamp() + 1
-            #time_arr[0, 4] = time_arr[0, 3] * time_float ** 4
+        if startTime_stamp > endTime_stamp:
-            #result_arr=result_arr_
+            raise
-            return result_arr_ # 位置矩阵           
+        delta_t = (endTime_stamp - startTime_stamp) / orbitnum
+        extractOrbits = []
+        #
+        temptime = startTime_stamp
+        while temptime < endTime_stamp:
+            temptime = temptime + delta_t
+            newOrbit = self.getTimeOrbitStamp(temptime)
+            extractOrbits.append(newOrbit)
+        newOrbit = self.getTimeOrbitStamp(endTime_stamp)
+        extractOrbits.append(newOrbit)
+        return extractOrbits  # 扩展的轨道节点
 
 
-def ReconstructionSatelliteOrbit(GPSPoints_list, starttime):
+def ReconstructionSatelliteOrbit(stateVectors, starttime):
     '''
     构建卫星轨道
     args:
         GPSPoints_list:卫星轨道点
         starttime:起算时间
     '''
-
+    GPSPoint_list = []
+    for sv in stateVectors:
+        GPSPoint = [sv.timeStamp.timestamp(),
+                    sv.xPosition,
+                    sv.yPosition,
+                    sv.zPosition,
+                    sv.xVelocity,
+                    sv.yVelocity,
+                    sv.zVelocity
+                    ]
+        GPSPoint_list.append(GPSPoint)
     SatelliteOrbitModel = SatelliteOrbitFitPoly()
-    if SatelliteOrbitModel.ReconstructionSatelliteOrbit(GPSPoints_list, starttime=starttime) is None:
+    if SatelliteOrbitModel.ReconstructionSatelliteOrbit(GPSPoint_list, starttime=starttime) is None:
         return None
-    return SatelliteOrbitModel
 
+    print("orbit test")
+    distance = []
+    for gpsPoint in GPSPoint_list:
+        temp_sv = SatelliteOrbitModel.SatelliteSpaceState(gpsPoint[0])
+        sv = np.array(temp_sv[1])
+        temp_distance = sv - gpsPoint[1:]
+        distance.append(temp_distance)
+    distance = np.array(distance)
+    print("orbit max:", np.max(distance))
+    print("orbit min:", np.min(distance))
+    return SatelliteOrbitModel
 
 
 ########
@@ -752,37 +914,40 @@ class GF3_SLC(Sensor):
         stopTime_s = dataStopTime.timestamp()
         centerTime_s = dataCenterTime.timestamp()
 
-        num = 0
+        tempOrbit = ReconstructionSatelliteOrbit(stateVectors, startTime_s)
-        time_list= []
+        As_arr = tempOrbit.A_arr.tolist()
-        for i in range(len(stateVectors)):
+        """ 构建轨道 """
-            # print(i,stateVectors[i].timeStamp)
+        self.frame.getOrbit().setsessionMode(self.product.satellite)
-            # tempOrbit.addvector(stateVectors[i].timeStamp.strftime("%Y-%m-%dT%H:%M:%S.%f"),
+        print("卫星型号")
-            #                     stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity,
+        print(self.frame.orbit.sessionMode)
-            #                     stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition)
+        self.frame.getOrbit().setPolyParams(tempOrbit.polynum, tempOrbit.starttime, As_arr)
-            if stateVectors[i].timeStamp>dataStartTime and stateVectors[i].timeStamp < dataStopTime:
+        self.frame.getOrbit().setDoppler(self.product.processInfo.DopplerParametersReferenceTime,
-                # tempOrbit.addvector(stateVectors[i].timeStamp.strftime("%Y-%m-%dT%H:%M:%S.%f"),
+                                         self.product.processInfo.DopplerRateValuesCoefficients)
-                #                 stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity,
+        """ """
-                #                 stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition)
+        newOrb = self.frame.getOrbit().getTimeOrbits(self.frame.sensingStart.strftime("%Y-%m-%dT%H:%M:%S.%f"),
-                # time_list.append(stateVectors[i].timeStamp)
+                                                     self.frame.sensingStop.strftime("%Y-%m-%dT%H:%M:%S.%f"),
-                time_list.append([stateVectors[i].timeStamp.timestamp(),
+                                                     orbitnum=1000)
-                          stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition,
+        for svect in newOrb:
-                          stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity])
+            self.frame.getOrbit().addStateVector(svect)
-                num += 1
+
-                # sv= StateVector()
+        # num = 0
-                # sv.setTime(stateVectors[i].timeStamp)
+        # time_list= []
-                # sv.setPosition([stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition])
+        # print("插值所有轨道")
-                # sv.setVelocity([stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity])
+        # for i in range(len(stateVectors)):
-                # self.frame.getOrbit().addStateVector(sv) # 插值结果写进了轨道模型中      
+        #
+        #     time_list.append([stateVectors[i].timeStamp.timestamp(),
+        #               stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition,
+        #               stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity])
         #     num += 1
-
+        #
-        model = ReconstructionSatelliteOrbit(time_list, starttime=centerTime_s)
+        # model = ReconstructionSatelliteOrbit(time_list, starttime=centerTime_s)
-        time_dif = ((stopTime_s + 30) - (startTime_s - 30)) / 1000
+        # time_dif = ((stopTime_s + 30) - (startTime_s - 30)) / 1000
-        time = np.zeros((1000, 1))
+        # time = np.zeros((1000, 1))
-        for i in range(1000):
+        # for i in range(1000):
-            time[i,:]=((startTime_s - 10) + time_dif * i)
+        #     time[i,:]=((startTime_s - 30) + time_dif * i)
-        t = time.reshape(-1)
+        # t = time.reshape(-1)
-    
+        #
-        statepoints = model.getSatelliteSpaceState(t)
+        # statepoints = model.getSatelliteSpaceState(t)
 
         # print("初始插值-----------------------------------------------------")
         # self.frame.setSensingStart(datetime.datetime.fromtimestamp(t[2]))
@@ -795,30 +960,14 @@ class GF3_SLC(Sensor):
         sensingMid = self.frame.sensingStart + datetime.timedelta(microseconds=int(diffTime*1e6))
         self.frame.setSensingMid(sensingMid)
         
-        # planet = self.frame.instrument.platform.planet
+        # for i, value in zip(range(len(statepoints)), statepoints):
-        # orbExt = OrbitExtender(planet=planet)
-        # orbExt.configure()
-        # newOrb = orbExt.extendOrbit(tempOrbit)
-        # tempOrbit.createOrbit() # 构建轨道模型
-
-        # newOrb=tempOrbit.getTimeOrbits(self.frame.sensingStart.strftime("%Y-%m-%dT%H:%M:%S.%f"), 
-        #                                 self.frame.sensingStop.strftime("%Y-%m-%dT%H:%M:%S.%f"),
-        #                                 orbitnum=500)
-        # for svect in newOrb:
         #     sv= StateVector()
-        #     sv.setTime(svect.UTCTime)
+        #     sv.setTime(datetime.datetime.fromtimestamp(t[i]))
-        #     sv.setPosition([svect.px,svect.py,svect.pz])
+        #     sv.setPosition([value[0],value[1],value[2]])
-        #     sv.setVelocity([svect.vx,svect.vy,svect.vz])
+        #     sv.setVelocity([value[3],value[4],value[5]])
         #     self.frame.getOrbit().addStateVector(sv) # 插值结果写进了轨道模型中
-        
+        #     # print("插值后的gps点", datetime.datetime.fromtimestamp(t[i]),value[0],value[1],value[2],value[3],value[4],value[5])
-        for i, value in zip(range(len(statepoints)), statepoints):
+        # print('Orbits list len is %d' %num)
-            sv= StateVector()
-            sv.setTime(datetime.datetime.fromtimestamp(t[i]))
-            sv.setPosition([value[0],value[1],value[2]])
-            sv.setVelocity([value[3],value[4],value[5]])
-            self.frame.getOrbit().addStateVector(sv) # 插值结果写进了轨道模型中
-            # print("插值后的gps点", datetime.datetime.fromtimestamp(t[i]),value[0],value[1],value[2],value[3],value[4],value[5])
-        print('Orbits list len is %d' %num)
         print('Successfully read state vectors from product XML')
 
     def extractPreciseOrbit(self, orbitfile, tstart, tend):

soilSalinity/SoilSalinity.xml

@@ -38,7 +38,7 @@
         <ParaType>File</ParaType>
         <DataType>tar.gz</DataType>
 		<ParaSource>Man</ParaSource>
-		<ParaValue>E:\202306hb\sar_img\GF3B_SYC_QPSI_008316_E116.2_N43.7_20230622_L1A_AHV_L10000202891-ortho.tar.gz</ParaValue>
+		<ParaValue>E:\2023xibei\GF3B_KSC_QPSI_010364_E86.0_N44.9_20231112_L1A_AHV_L10000263358-ortho.tar.gz</ParaValue>
         <MaxValue>DEFAULT</MaxValue>
         <MinValue>DEFAULT</MinValue>
         <OptionValue>DEFAULT</OptionValue>
@@ -66,7 +66,7 @@
         <ParaType>File</ParaType>
         <DataType>csv</DataType>
 		<ParaSource>Man</ParaSource>
-        <ParaValue>E:\202306pj\result2\soilSality(1).csv</ParaValue>
+        <ParaValue>E:\2023xibei\soilSality.csv</ParaValue>
         <MaxValue>DEFAULT</MaxValue>
         <MinValue>DEFAULT</MinValue>
         <OptionValue>DEFAULT</OptionValue>
@@ -138,7 +138,7 @@
         <ParaType>File</ParaType>
         <DataType>tar.gz</DataType>
 		<ParaSource>Man</ParaSource>
-        <ParaValue>D:\micro\WorkSpace\SoilSalinity\Output\GF3B_SYC_QPSI_008316_E116.2_N43.7_20230622_L1A_AHV_L10000202891-ortho-SSAA.tar.gz</ParaValue>
+        <ParaValue>D:\micro\WorkSpace\SoilSalinity\Output\GF3B_KSC_QPSI_010364_E86.0_N44.9_20231112_L1A_AHV_L10000263358-ortho-SSAA.tar.gz</ParaValue>
         <MaxValue>DEFAULT</MaxValue>
         <MinValue>DEFAULT</MinValue>
         <OptionValue>DEFAULT</OptionValue>

soilSalinity/SoilSalinityMain.py

@@ -101,7 +101,7 @@ class SalinityMain:
         self.__workspace_path = self.__alg_xml_handler.get_workspace_path()
         self.__create_work_space()
 
-        self.__processing_paras = InitPara.init_processing_paras(self.__input_paras)
+        self.__processing_paras = InitPara.init_processing_paras(self.__input_paras, self.__workspace_preprocessed_path)
         self.__processing_paras.update(self.get_tar_gz_inf(self.__processing_paras["sar_path0"]))
 
         SrcImageName = os.path.split(self.__input_paras["AHV"]['ParaValue'])[1].split('.tar.gz')[0]
@@ -159,7 +159,7 @@ class SalinityMain:
         """
         预处理
         """
-        para_names_geo = ["Covering", "NDVI", 'sim_ori']
+        para_names_geo = ["Covering", "NDVI", 'sim_ori']  #
         p = pp()
         p.check_img_projection(self.__workspace_preprocessing_path, para_names_geo, self.__processing_paras)
         #计算roi
@@ -360,15 +360,15 @@ class SalinityMain:
         block_features_dir, block_features_name = bp.get_file_names(self.__workspace_block_tif_processed_path + 'features\\', ['tif'])
         X_train, Y_train = ml.gene_train_data(block_features_dir, rows, cols, block_size, measured_data_img)
 
-        optimal_feature = ml.sel_optimal_feature_set(X_train, Y_train, threshold=0.01)
+        # optimal_feature = ml.sel_optimal_feature_set(X_train, Y_train, threshold=0.01)
-        optimal_feature = ml.remove_correlation_feature(X_train, optimal_feature, threshold=0.85)
+        # optimal_feature = ml.remove_correlation_feature(X_train, optimal_feature, threshold=0.85)
-        X_train = X_train[:, optimal_feature]
+        # X_train = X_train[:, optimal_feature]
 
         feature_name_list = []
         for name in dir_dict.keys():
             feature_name_list.append(name)
         logger.info('feature_list:%s', feature_name_list)
-        logger.info('train_feature:%s', np.array(feature_name_list)[optimal_feature])
+        # logger.info('train_feature:%s', np.array(feature_name_list)[optimal_feature])
         logger.info('generating training set success!')
         logger.info('progress bar: 80%')
 
@@ -383,7 +383,7 @@ class SalinityMain:
         for path, name, n in zip(block_features_dir, block_features_name, range(len(block_features_dir))):
             features_array = self.imageHandler.get_data(path)
             X_test = np.reshape(features_array, (features_array.shape[0], features_array[0].size)).T
-            X_test = X_test[:, optimal_feature]
+            # X_test = X_test[:, optimal_feature]
             Y_test = pls.predict(X_test)
             Y_test[Y_test < soil_salinity_value_min] = soil_salinity_value_min
             Y_test[Y_test > soil_salinity_value_max] = soil_salinity_value_max

tool/algorithm/ml/machineLearning.py

@@ -373,16 +373,19 @@ class MachineLeaning:
             block_col = col//block_size
 
             if block_row == block_rows-1:
-                part_img_row = row - (rows - block_size)
+                # part_img_row = row - (rows - block_size)
+                part_img_row = row - (block_row * block_size)
             else:
                 part_img_row = row % block_size
 
             if block_col == block_cols-1:
-                part_img_col = col - (cols-block_size)
+                # part_img_col = col - (cols-block_size)
+                part_img_col = col - (block_col * block_size)
             else:
                 part_img_col = col % block_size
 
-            features_path = block_features_dir[block_row*(block_rows-1) + block_col]
+            # features_path = block_features_dir[block_row*(block_rows-1) + block_col]
+            features_path = block_features_dir[block_row*block_cols + block_col]
             features_array = ImageHandler().get_data(features_path)
 
             feature = features_array[:, part_img_row, part_img_col]

tool/algorithm/transforml1a/transHandle.py

@@ -215,14 +215,14 @@ class TransImgL1A:
                     c_count = c_max - c_min  # 列数
 
                     #
-                    mask_l1a = np.zeros((r_count, c_count)) * np.nan  # 创建目标大小的行列号
+                    mask_l1a = np.zeros((int(r_count), int(c_count))) * np.nan  # 创建目标大小的行列号
                     mask = SAR_GEO.gereratorMask(r_idx.astype(np.float64), c_idx.astype(np.float64).astype(np.float64),
                                                  mask_l1a)  # 这个函数修改了
 
-                    self._begin_r = r_min
+                    self._begin_r = int(r_min)
-                    self._end_r = r_max
+                    self._end_r = int(r_max)
-                    self._begin_c = c_min
+                    self._begin_c = int(c_min)
-                    self._end_c = c_max
+                    self._end_c = int(c_max)
                     self._mask = mask
 
     def cut_L1A(self, in_path, out_path):