diff --git a/Ortho/Ortho.xml b/Ortho/Ortho.xml index 266cd518..47738901 100644 --- a/Ortho/Ortho.xml +++ b/Ortho/Ortho.xml @@ -1,7 +1,7 @@ CSAR_202107275419_0001-0 - E:\Result_GF3\ + D:\micro\WorkSpace\ File ElementAlg @@ -45,7 +45,7 @@ File tar.gz Cal - E:\GF3Data\soilMoisture\GF3B_SYC_QPSI_008316_E116.1_N43.3_20230622_L1A_AHV_L10000202892.tar.gz + D:\micro\microproduct_depdence\GF3-Deformation\download\cls\GF3_SAY_FSI_001614_E113.2_N34.5_20161129_L1A_HHHV_L10002015686.tar.gz True False File @@ -58,9 +58,9 @@ DEM数字高程影像 30m分辨率DEM数字高程影像tif File - zip + File Cal - 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 + D:\micro\microproduct_depdence\GF3-Deformation\dem True True File @@ -92,7 +92,7 @@ File tar.gz Cal - E:\Result_GF3\Ortho\Output\GF3B_SYC_QPSI_008316_E116.1_N43.3_20230622_L1A_AHV_L10000202892-ortho.tar.gz + D:\micro\WorkSpace\Ortho\Output\GF3_SAY_FSI_001614_E113.2_N34.5_20161129_L1A_HHHV_L10002015686-ortho.tar.gz DEFAULT DEFAULT DEFAULT diff --git a/dem-C-SAR/Dem.xml b/dem-C-SAR/Dem.xml index 66d1210d..9f8852df 100644 --- a/dem-C-SAR/Dem.xml +++ b/dem-C-SAR/Dem.xml @@ -71,7 +71,7 @@ value string Man - 34.60;34.67;113.05;113.18 + 34.6;34.63;113.1;113.15 True True UploadInput diff --git a/dem-C-SAR/DemMain.py b/dem-C-SAR/DemMain.py index 8dbe285d..95c8026f 100644 --- a/dem-C-SAR/DemMain.py +++ b/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)) diff --git a/dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Orbit/Orbit.py b/dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Orbit/Orbit.py index dcaafd15..e055cfa2 100644 --- a/dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Orbit/Orbit.py +++ b/dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Orbit/Orbit.py @@ -234,7 +234,16 @@ class Orbit(Component): self._orbitSource = source or None self._referenceFrame = None self._stateVectors.configure() - #self._stateVectors = stateVectors or [] + """ 自定义 """ + 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) return None @@ -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,38 +448,46 @@ class Orbit(Component): @raises NotImplementedError: if the desired interpolation method cannot be decoded """ - if time not in self: - raise ValueError( - "Time stamp (%s) falls outside of the interpolation interval [%s:%s]" % - (time,self.minTime,self.maxTime) + 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) ) - if method == 'linear': - newSV = self._linearOrbitInterpolation(time) - elif method == 'legendre': - newSV = self._legendreOrbitInterpolation(time) - elif method == 'hermite': - newSV = self._hermiteOrbitInterpolation(time) - else: - raise NotImplementedError( - "Orbit interpolation type %s, is not implemented" % method + if method == 'linear': + newSV = self._linearOrbitInterpolation(time) + elif method == 'legendre': + newSV = self._legendreOrbitInterpolation(time) + elif method == 'hermite': + newSV = self._hermiteOrbitInterpolation(time) + else: + raise NotImplementedError( + "Orbit interpolation type %s, is not implemented" % method ) - return newSV + 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 time not in self: - raise ValueError("Time stamp (%s) falls outside of the interpolation interval [%s:%s]" - % (time,self.minTime,self.maxTime)) - try: - return getattr(self, '_'+method+'OrbitInterpolation')(time) - except AttributeError: - pass - raise NotImplementedError( - "Orbit interpolation type %s, is not implemented" % method + 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)) + try: + return getattr(self, '_' + method + 'OrbitInterpolation')(time) + except AttributeError: + pass + 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. @@ -1011,51 +1091,183 @@ class Orbit(Component): wvl = 0.0 if planet is None: - refElp = Planet(pname='Earth'). ellipsoid + 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) - 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 + 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) - pos = np.array(sv.getPosition()) - vel = np.array(sv.getVelocity()) + 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 - dr = xyz-pos - rng = np.linalg.norm(dr) + pos1 = np.array(sv.getPosition()) # 卫星坐标 + vel1 = np.array(sv.getVelocity()) # 卫星速度 + dr1 = pos1 - xyz + rng1 = np.linalg.norm(dr1) # 斜距 - 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 + # ((R_s1.array() * V_s1.array()).rowwise().sum().array() * (-2) / (R * this->lamda))[0]; + FdTheory1 = -2 / (rng1 * wvl) * np.dot(dr1, vel1) - fn = dopfact - fdop * rng - c1 = -np.dot(vel, vel) - c2 = (fdop/rng + fdopder) + 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 - fnprime = c1 + c2 * dopfact + 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 - tguess = tguess - datetime.timedelta(seconds = fn/fnprime) - if abs(tguess - t_prev_guess).total_seconds() < 5e-9: - break - else: - t_prev_guess = tguess + 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 outOfBounds: - raise Exception('Interpolation time out of bounds') + 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 + for ii in range(51): + try: + sv = self.interpolateOrbit(tguess, method='hermite') + 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) # 求斜距 + # print("rng", rng) + + dopfact = np.dot(dr, vel) # fd 公式 + # print("dopfact", dopfact) + + fdop = doppler(DTU.seconds_since_midnight(tguess), rng) * wvl * 0.5 + # 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) + # print("输出的tguess", tguess) + # print(outOfBounds) + print("end ------------------------------------------------------------\n") + 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 diff --git a/dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Sensor/GF3_SLC.py b/dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Sensor/GF3_SLC.py index 54a4b44c..3c48a8ea 100644 --- a/dem-C-SAR/ISCEApp/_internal/isce/components/isceobj/Sensor/GF3_SLC.py +++ b/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类，为拟合多项式实现方法 @@ -171,75 +325,58 @@ class SatelliteOrbitFitPoly(SatelliteOrbit): def __init__(self) -> None: super().__init__() - self.modelName="多项式" - self.polynum=4 + self.modelName = "多项式" + self.polynum = 4 def ReconstructionSatelliteOrbit(self, GPSPoints_list, starttime): - if len(GPSPoints_list)==2: - self.polynum=1 + 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) # 记录时间坐标 + 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=[] + 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()) + 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) # 记录时间坐标 + A_arr = np.zeros((self.polynum, 6), 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]) - state_arr[i, :] = np.array(GPSPoint[1:],dtype=np.float64).reshape(1,6) # 空间坐标 - self.model_f=[] + time_ = GPSPoint[0] - self.starttime # 为了保证精度，对时间进行缩放 + 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()) - ''' 测试误差 - 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() - ''' + 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 = A_arr.copy() return self.A_arr else: self.A_arr = None @@ -257,66 +394,91 @@ class SatelliteOrbitFitPoly(SatelliteOrbit): if self.model_f is None: return None - result_arr=np.zeros((1,7)) + 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_ # 位置矩阵 + # + 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 + + return [time_float, [px, py, pz, vx, vy, vz]] + + def getTimeOrbitStamp(self, UTCStartTime_float): + sv = _StateVector() + temp_sv = self.SatelliteSpaceState(UTCStartTime_float) + sv.timeStamp = datetime.datetime.fromtimestamp(UTCStartTime_float) + sv.xPosition = temp_sv[1][0, 0] + sv.yPosition = temp_sv[1][0, 1] + sv.zPosition = temp_sv[1][0, 2] + sv.xVelocity = temp_sv[1][0, 3] + sv.yVelocity = temp_sv[1][0, 4] + sv.zVelocity = temp_sv[1][0, 5] + return sv + + def getTimeOrbits(self, UTCStartTime, UTCEndTime, orbitnum=1000): + # + startTime_stamp = datetime.datetime.strptime(UTCStartTime, "%Y-%m-%dT%H:%M:%S.%f").timestamp() - 1 + endTime_stamp = datetime.datetime.strptime(UTCEndTime, "%Y-%m-%dT%H:%M:%S.%f").timestamp() + 1 + if startTime_stamp > endTime_stamp: + raise + 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 ######## @@ -621,7 +783,7 @@ class GF3_SLC(Sensor): print("effective PRF %f" % (prf) ) print('-----------------------------------------------------') - lineFlip = (self.product.direction.upper() == 'DECREASING') + lineFlip = (self.product.direction.upper() == 'DECREASING') dataStartTime = self.product.imageInfo.imagingStartTime dataStopTime = self.product.imageInfo.imagingEndTime @@ -750,39 +912,42 @@ class GF3_SLC(Sensor): startTime_s = dataStartTime.timestamp() stopTime_s = dataStopTime.timestamp() - centerTime_s = dataCenterTime.timestamp() + centerTime_s = dataCenterTime.timestamp() - num = 0 - time_list= [] - for i in range(len(stateVectors)): - # print(i,stateVectors[i].timeStamp) - # tempOrbit.addvector(stateVectors[i].timeStamp.strftime("%Y-%m-%dT%H:%M:%S.%f"), - # stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity, - # stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition) - if stateVectors[i].timeStamp>dataStartTime and stateVectors[i].timeStamp < dataStopTime: - # tempOrbit.addvector(stateVectors[i].timeStamp.strftime("%Y-%m-%dT%H:%M:%S.%f"), - # stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity, - # stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition) - # time_list.append(stateVectors[i].timeStamp) - 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 - # sv= StateVector() - # sv.setTime(stateVectors[i].timeStamp) - # sv.setPosition([stateVectors[i].xPosition, stateVectors[i].yPosition, stateVectors[i].zPosition]) - # sv.setVelocity([stateVectors[i].xVelocity, stateVectors[i].yVelocity, stateVectors[i].zVelocity]) - # self.frame.getOrbit().addStateVector(sv) # 插值结果写进了轨道模型中 - # num+=1 + tempOrbit = ReconstructionSatelliteOrbit(stateVectors, startTime_s) + As_arr = tempOrbit.A_arr.tolist() + """ 构建轨道 """ + self.frame.getOrbit().setsessionMode(self.product.satellite) + print("卫星型号") + print(self.frame.orbit.sessionMode) + self.frame.getOrbit().setPolyParams(tempOrbit.polynum, tempOrbit.starttime, As_arr) + self.frame.getOrbit().setDoppler(self.product.processInfo.DopplerParametersReferenceTime, + self.product.processInfo.DopplerRateValuesCoefficients) + """ """ + newOrb = self.frame.getOrbit().getTimeOrbits(self.frame.sensingStart.strftime("%Y-%m-%dT%H:%M:%S.%f"), + self.frame.sensingStop.strftime("%Y-%m-%dT%H:%M:%S.%f"), + orbitnum=1000) + for svect in newOrb: + self.frame.getOrbit().addStateVector(svect) - model = ReconstructionSatelliteOrbit(time_list, starttime=centerTime_s) - time_dif = ((stopTime_s + 30) - (startTime_s - 30)) / 1000 - time = np.zeros((1000, 1)) - for i in range(1000): - time[i,:]=((startTime_s - 10) + time_dif * i) - t = time.reshape(-1) - - statepoints = model.getSatelliteSpaceState(t) + # num = 0 + # time_list= [] + # print("插值所有轨道") + # for i in range(len(stateVectors)): + # + # 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) + # time_dif = ((stopTime_s + 30) - (startTime_s - 30)) / 1000 + # time = np.zeros((1000, 1)) + # for i in range(1000): + # time[i,:]=((startTime_s - 30) + time_dif * i) + # t = time.reshape(-1) + # + # statepoints = model.getSatelliteSpaceState(t) # print("初始插值-----------------------------------------------------") # self.frame.setSensingStart(datetime.datetime.fromtimestamp(t[2])) @@ -794,31 +959,15 @@ class GF3_SLC(Sensor): diffTime = DTUtil.timeDeltaToSeconds(self.frame.sensingStop-self.frame.sensingStart)/2.0 sensingMid = self.frame.sensingStart + datetime.timedelta(microseconds=int(diffTime*1e6)) self.frame.setSensingMid(sensingMid) - - # planet = self.frame.instrument.platform.planet - # 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.setPosition([svect.px,svect.py,svect.pz]) - # sv.setVelocity([svect.vx,svect.vy,svect.vz]) - # self.frame.getOrbit().addStateVector(sv) # 插值结果写进了轨道模型中 - for i, value in zip(range(len(statepoints)), statepoints): - 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) + # for i, value in zip(range(len(statepoints)), statepoints): + # 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): diff --git a/dem-C-SAR/ISCEApp/stackStripMap.exe b/dem-C-SAR/ISCEApp/stackStripMap.exe index 07f287d8..ac869922 100644 Binary files a/dem-C-SAR/ISCEApp/stackStripMap.exe and b/dem-C-SAR/ISCEApp/stackStripMap.exe differ diff --git a/dem-C-SAR/ISCEApp/stripmapWrapper.exe b/dem-C-SAR/ISCEApp/stripmapWrapper.exe index 0f165456..46b72fb3 100644 Binary files a/dem-C-SAR/ISCEApp/stripmapWrapper.exe and b/dem-C-SAR/ISCEApp/stripmapWrapper.exe differ diff --git a/soilSalinity/SoilSalinity.xml b/soilSalinity/SoilSalinity.xml index 3bb6fd9a..a5b51787 100644 --- a/soilSalinity/SoilSalinity.xml +++ b/soilSalinity/SoilSalinity.xml @@ -38,7 +38,7 @@ File tar.gz Man - E:\202306hb\sar_img\GF3B_SYC_QPSI_008316_E116.2_N43.7_20230622_L1A_AHV_L10000202891-ortho.tar.gz + E:\2023xibei\GF3B_KSC_QPSI_010364_E86.0_N44.9_20231112_L1A_AHV_L10000263358-ortho.tar.gz DEFAULT DEFAULT DEFAULT @@ -66,7 +66,7 @@ File csv Man - E:\202306pj\result2\soilSality(1).csv + E:\2023xibei\soilSality.csv DEFAULT DEFAULT DEFAULT @@ -138,7 +138,7 @@ File tar.gz Man - D:\micro\WorkSpace\SoilSalinity\Output\GF3B_SYC_QPSI_008316_E116.2_N43.7_20230622_L1A_AHV_L10000202891-ortho-SSAA.tar.gz + D:\micro\WorkSpace\SoilSalinity\Output\GF3B_KSC_QPSI_010364_E86.0_N44.9_20231112_L1A_AHV_L10000263358-ortho-SSAA.tar.gz DEFAULT DEFAULT DEFAULT diff --git a/soilSalinity/SoilSalinityMain.py b/soilSalinity/SoilSalinityMain.py index 2090b787..eeb948f3 100644 --- a/soilSalinity/SoilSalinityMain.py +++ b/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.remove_correlation_feature(X_train, optimal_feature, threshold=0.85) - X_train = X_train[:, optimal_feature] + # 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) + # 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 diff --git a/tool/algorithm/ml/machineLearning.py b/tool/algorithm/ml/machineLearning.py index 06abfcc3..eac2e261 100644 --- a/tool/algorithm/ml/machineLearning.py +++ b/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] diff --git a/tool/algorithm/transforml1a/transHandle.py b/tool/algorithm/transforml1a/transHandle.py index f53ce94f..04a7a30c 100644 --- a/tool/algorithm/transforml1a/transHandle.py +++ b/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._end_r = r_max - self._begin_c = c_min - self._end_c = c_max + self._begin_r = int(r_min) + self._end_r = int(r_max) + self._begin_c = int(c_min) + self._end_c = int(c_max) self._mask = mask def cut_L1A(self, in_path, out_path):