From 3fa7abc3fac0c7e311980393362bcd9155337156 Mon Sep 17 00:00:00 2001 From: Yang Lei Date: Tue, 26 Oct 2021 08:18:18 -0700 Subject: [PATCH 1/2] fix bugs for calculating error and stable shift in nc packaging --- netcdf_output.py | 39 +++++++++++--------- testautoRIFT.py | 86 +++++++++++++++++++++++++++----------------- testautoRIFT_ISCE.py | 86 +++++++++++++++++++++++++++----------------- 3 files changed, 130 insertions(+), 81 deletions(-) diff --git a/netcdf_output.py b/netcdf_output.py index 17d0cad..e398c1a 100755 --- a/netcdf_output.py +++ b/netcdf_output.py @@ -227,7 +227,7 @@ def netCDF_read_intermediate(filename='./autoRIFT_intermediate.nc'): return Dx, Dy, InterpMask, ChipSizeX, GridSpacingX, ScaleChipSizeY, SearchLimitX, SearchLimitY, origSize, noDataMask -def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, offset2vr, offset2va, MM, VXref, VYref, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector): +def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, offset2vr, offset2va, MM, VXref, VYref, DXref, DYref, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector): vx_mean_shift = offset2vx_1 * dx_mean_shift + offset2vx_2 * dy_mean_shift @@ -266,6 +266,11 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 VR = VR.astype(np.float32) VA = VA.astype(np.float32) + VRref = DXref * offset2vr + VAref = (DYref) * offset2va + VRref = VRref.astype(np.float32) + VAref = VAref.astype(np.float32) + # vr_mean_shift = dx_mean_shift * rangePixelSize / dt * 365.0 * 24.0 * 3600.0 # va_mean_shift = (dy_mean_shift) * azimuthPixelSize / dt * 365.0 * 24.0 * 3600.0 @@ -620,14 +625,14 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 if stable_count != 0: - temp = VX.copy() + temp = VX.copy() - VXref.copy() temp[np.logical_not(SSM)] = np.nan # vx_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) vx_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) else: vx_error_mask = np.nan if stable_count1 != 0: - temp = VX.copy() + temp = VX.copy() - VXref.copy() temp[np.logical_not(SSM1)] = np.nan # vx_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) vx_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) @@ -703,14 +708,14 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 if stable_count != 0: - temp = VY.copy() + temp = VY.copy() - VYref.copy() temp[np.logical_not(SSM)] = np.nan # vy_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) vy_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) else: vy_error_mask = np.nan if stable_count1 != 0: - temp = VY.copy() + temp = VY.copy() - VYref.copy() temp[np.logical_not(SSM1)] = np.nan # vy_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) vy_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) @@ -834,14 +839,14 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 var.setncattr('units','m/y') if stable_count != 0: - temp = VR.copy() + temp = VR.copy() - VRref.copy() temp[np.logical_not(SSM)] = np.nan # vr_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) vr_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) else: vr_error_mask = np.nan if stable_count1 != 0: - temp = VR.copy() + temp = VR.copy() - VRref.copy() temp[np.logical_not(SSM1)] = np.nan # vr_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) vr_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) @@ -914,14 +919,14 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 if stable_count != 0: - temp = VA.copy() + temp = VA.copy() - VAref.copy() temp[np.logical_not(SSM)] = np.nan # va_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) va_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) else: va_error_mask = np.nan if stable_count1 != 0: - temp = VA.copy() + temp = VA.copy() - VAref.copy() temp[np.logical_not(SSM1)] = np.nan # va_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) va_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) @@ -985,22 +990,22 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 # fuse the (slope parallel & reference) flow-based range-projected result with the raw observed range/azimuth-based result if stable_count_p != 0: - temp = VXP.copy() + temp = VXP.copy() - VXref.copy() temp[np.logical_not(SSM)] = np.nan # vxp_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) vxp_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) - temp = VYP.copy() + temp = VYP.copy() - VYref.copy() temp[np.logical_not(SSM)] = np.nan # vyp_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) vyp_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) if stable_count1_p != 0: - temp = VXP.copy() + temp = VXP.copy() - VXref.copy() temp[np.logical_not(SSM1)] = np.nan # vxp_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) vxp_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) - temp = VYP.copy() + temp = VYP.copy() - VYref.copy() temp[np.logical_not(SSM1)] = np.nan # vyp_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) vyp_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) @@ -1088,14 +1093,14 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 if stable_count_p != 0: - temp = VXP.copy() + temp = VXP.copy() - VXref.copy() temp[np.logical_not(SSM)] = np.nan # vxp_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) vxp_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) else: vxp_error_mask = np.nan if stable_count1_p != 0: - temp = VXP.copy() + temp = VXP.copy() - VXref.copy() temp[np.logical_not(SSM1)] = np.nan # vxp_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) vxp_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) @@ -1170,14 +1175,14 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 if stable_count_p != 0: - temp = VYP.copy() + temp = VYP.copy() - VYref.copy() temp[np.logical_not(SSM)] = np.nan # vyp_error_mask = np.std(temp[(temp > -500)&(temp < 500)]) vyp_error_mask = np.std(temp[np.logical_not(np.isnan(temp))]) else: vyp_error_mask = np.nan if stable_count1_p != 0: - temp = VYP.copy() + temp = VYP.copy() - VYref.copy() temp[np.logical_not(SSM1)] = np.nan # vyp_error_slow = np.std(temp[(temp > -500)&(temp < 500)]) vyp_error_slow = np.std(temp[np.logical_not(np.isnan(temp))]) diff --git a/testautoRIFT.py b/testautoRIFT.py index 319bf2b..aaf3c94 100755 --- a/testautoRIFT.py +++ b/testautoRIFT.py @@ -707,7 +707,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # stable_count = np.sum(SSM & np.logical_not(np.isnan(DX)) & (DX-DXref > -5) & (DX-DXref < 5) & (DY-DYref > -5) & (DY-DYref < 5)) stable_count = np.sum(SSM & np.logical_not(np.isnan(DX))) - V_temp = np.sqrt(VX**2 + VY**2) + V_temp = np.sqrt(VXref**2 + VYref**2) try: V_temp_threshold = np.percentile(V_temp[np.logical_not(np.isnan(V_temp))],25) SSM1 = (V_temp <= V_temp_threshold) @@ -810,18 +810,21 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search - from datetime import datetime, timedelta -# d0 = datetime(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) -# d1 = datetime(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) + from datetime import date, datetime +# d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) +# d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) d0 = datetime.strptime(master_dt,"%Y%m%dT%H:%M:%S.%f") d1 = datetime.strptime(slave_dt,"%Y%m%dT%H:%M:%S.%f") - date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() - date_dt = np.float64(date_dt_base) + date_dt_base = d1 - d0 + date_dt = np.float64(date_dt_base.days) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'time_standard_img1':'UTC','absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'flight_direction_img1':flight_direction_m,'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'time_standard_img2':'UTC','absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'flight_direction_img2':flight_direction_s,'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([[0.0356, 0.0501, 0.0266, 0.0622, 0.0357, 0.0501], @@ -830,7 +833,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search netcdf_file = no.netCDF_packaging( VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, offset2vr, offset2va, MM, VXref, VYref, - rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, + DXref, DYref, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector ) @@ -863,6 +866,12 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # # master_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+master_MTL_path))[2][1:-2],':') # slave_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+slave_MTL_path))[2][1:-2],':') + master_time = ['00','00','00'] + slave_time = ['00','00','00'] + + from datetime import time as time1 + master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) + slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) import netcdf_output as no pair_type = 'optical' @@ -882,19 +891,22 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import datetime, timedelta - d0 = datetime(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) - d1 = datetime(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) - date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() - date_dt = np.float64(date_dt_base) + from datetime import date + d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) + d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(date_dt_base.days) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%d") - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - - master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + master_dt = master_split[3][0:8] + master_time.strftime("T%H:%M:%S") + slave_dt = slave_split[3][0:8] + slave_time.strftime("T%H:%M:%S") IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} @@ -903,7 +915,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search netcdf_file = no.netCDF_packaging( VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, None, None, MM, VXref, VYref, - XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, + None, None, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector ) @@ -931,6 +943,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_filename = master_split[0][-3:]+'_'+master_split[2]+'_'+master_split[4][:3]+'_'+os.path.basename(master_path) slave_filename = slave_split[0][-3:]+'_'+slave_split[2]+'_'+slave_split[4][:3]+'_'+os.path.basename(slave_path) + master_time = ['00','00','00'] + slave_time = ['00','00','00'] + + from datetime import time as time1 + master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) + slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) + import netcdf_output as no pair_type = 'optical' detection_method = 'feature' @@ -947,19 +966,22 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import datetime, timedelta - d0 = datetime(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) - d1 = datetime(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) - date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() - date_dt = np.float64(date_dt_base) + from datetime import date + d0 = date(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) + d1 = date(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(date_dt_base.days) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%d") - date_ct = d0 + (d1 - d0) / 2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - - master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + master_dt = master_split[2] + master_time.strftime("T%H:%M:%S") + slave_dt = slave_split[2] + slave_time.strftime("T%H:%M:%S") IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} @@ -968,7 +990,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search netcdf_file = no.netCDF_packaging( VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, None, None, MM, VXref, VYref, - XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, + None, None, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector ) diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py index 6f08b92..e98a6f1 100755 --- a/testautoRIFT_ISCE.py +++ b/testautoRIFT_ISCE.py @@ -712,7 +712,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # stable_count = np.sum(SSM & np.logical_not(np.isnan(DX)) & (DX-DXref > -5) & (DX-DXref < 5) & (DY-DYref > -5) & (DY-DYref < 5)) stable_count = np.sum(SSM & np.logical_not(np.isnan(DX))) - V_temp = np.sqrt(VX**2 + VY**2) + V_temp = np.sqrt(VXref**2 + VYref**2) try: V_temp_threshold = np.percentile(V_temp[np.logical_not(np.isnan(V_temp))],25) SSM1 = (V_temp <= V_temp_threshold) @@ -815,18 +815,21 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search - from datetime import datetime, timedelta -# d0 = datetime(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) -# d1 = datetime(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) + from datetime import date, datetime +# d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) +# d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) d0 = datetime.strptime(master_dt,"%Y%m%dT%H:%M:%S.%f") d1 = datetime.strptime(slave_dt,"%Y%m%dT%H:%M:%S.%f") - date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() - date_dt = np.float64(date_dt_base) + date_dt_base = d1 - d0 + date_dt = np.float64(date_dt_base.days) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'time_standard_img1':'UTC','absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'flight_direction_img1':flight_direction_m,'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'time_standard_img2':'UTC','absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'flight_direction_img2':flight_direction_s,'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([[0.0356, 0.0501, 0.0266, 0.0622, 0.0357, 0.0501], @@ -835,7 +838,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search netcdf_file = no.netCDF_packaging( VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, offset2vr, offset2va, MM, VXref, VYref, - rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, + DXref, DYref, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector ) @@ -868,6 +871,12 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # # master_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+master_MTL_path))[2][1:-2],':') # slave_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+slave_MTL_path))[2][1:-2],':') + master_time = ['00','00','00'] + slave_time = ['00','00','00'] + + from datetime import time as time1 + master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) + slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) import netcdf_output as no pair_type = 'optical' @@ -886,19 +895,22 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import datetime, timedelta - d0 = datetime(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) - d1 = datetime(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) - date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() - date_dt = np.float64(date_dt_base) + from datetime import date + d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) + d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(date_dt_base.days) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%d") - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - - master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + master_dt = master_split[3][0:8] + master_time.strftime("T%H:%M:%S") + slave_dt = slave_split[3][0:8] + slave_time.strftime("T%H:%M:%S") IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} @@ -907,7 +919,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search netcdf_file = no.netCDF_packaging( VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, None, None, MM, VXref, VYref, - XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, + None, None, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector ) @@ -935,6 +947,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_filename = master_split[0][-3:]+'_'+master_split[2]+'_'+master_split[4][:3]+'_'+os.path.basename(master_path) slave_filename = slave_split[0][-3:]+'_'+slave_split[2]+'_'+slave_split[4][:3]+'_'+os.path.basename(slave_path) + master_time = ['00','00','00'] + slave_time = ['00','00','00'] + + from datetime import time as time1 + master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) + slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) + import netcdf_output as no pair_type = 'optical' detection_method = 'feature' @@ -951,19 +970,22 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import datetime, timedelta - d0 = datetime(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) - d1 = datetime(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) - date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() - date_dt = np.float64(date_dt_base) + from datetime import date + d0 = date(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) + d1 = date(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(date_dt_base.days) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%d") - date_ct = d0 + (d1 - d0) / 2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - - master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + master_dt = master_split[2] + master_time.strftime("T%H:%M:%S") + slave_dt = slave_split[2] + slave_time.strftime("T%H:%M:%S") IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} @@ -972,7 +994,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search netcdf_file = no.netCDF_packaging( VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, None, None, MM, VXref, VYref, - XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, + None, None, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_count1, stable_shift_applied, dx_mean_shift, dy_mean_shift, dx_mean_shift1, dy_mean_shift1, error_vector ) From 12d37070fa4aeedd250355d463b13e5c178bd9e6 Mon Sep 17 00:00:00 2001 From: Yang Lei Date: Tue, 26 Oct 2021 08:26:24 -0700 Subject: [PATCH 2/2] fix bugs for calculating error and stable shift in nc packaging --- testautoRIFT.py | 78 ++++++++++++++++---------------------------- testautoRIFT_ISCE.py | 78 ++++++++++++++++---------------------------- 2 files changed, 56 insertions(+), 100 deletions(-) diff --git a/testautoRIFT.py b/testautoRIFT.py index aaf3c94..c9f8e0c 100755 --- a/testautoRIFT.py +++ b/testautoRIFT.py @@ -810,21 +810,18 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search - from datetime import date, datetime -# d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) -# d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) + from datetime import datetime, timedelta +# d0 = datetime(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) +# d1 = datetime(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) d0 = datetime.strptime(master_dt,"%Y%m%dT%H:%M:%S.%f") d1 = datetime.strptime(slave_dt,"%Y%m%dT%H:%M:%S.%f") - date_dt_base = d1 - d0 - date_dt = np.float64(date_dt_base.days) + date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() + date_dt = np.float64(date_dt_base) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - else: - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'time_standard_img1':'UTC','absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'flight_direction_img1':flight_direction_m,'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'time_standard_img2':'UTC','absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'flight_direction_img2':flight_direction_s,'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([[0.0356, 0.0501, 0.0266, 0.0622, 0.0357, 0.0501], @@ -866,12 +863,6 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # # master_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+master_MTL_path))[2][1:-2],':') # slave_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+slave_MTL_path))[2][1:-2],':') - master_time = ['00','00','00'] - slave_time = ['00','00','00'] - - from datetime import time as time1 - master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) - slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) import netcdf_output as no pair_type = 'optical' @@ -891,22 +882,19 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import date - d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) - d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(date_dt_base.days) + from datetime import datetime, timedelta + d0 = datetime(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) + d1 = datetime(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) + date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() + date_dt = np.float64(date_dt_base) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") - else: - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%d") - master_dt = master_split[3][0:8] + master_time.strftime("T%H:%M:%S") - slave_dt = slave_split[3][0:8] + slave_time.strftime("T%H:%M:%S") + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + + master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} @@ -943,13 +931,6 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_filename = master_split[0][-3:]+'_'+master_split[2]+'_'+master_split[4][:3]+'_'+os.path.basename(master_path) slave_filename = slave_split[0][-3:]+'_'+slave_split[2]+'_'+slave_split[4][:3]+'_'+os.path.basename(slave_path) - master_time = ['00','00','00'] - slave_time = ['00','00','00'] - - from datetime import time as time1 - master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) - slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) - import netcdf_output as no pair_type = 'optical' detection_method = 'feature' @@ -966,22 +947,19 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import date - d0 = date(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) - d1 = date(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(date_dt_base.days) + from datetime import datetime, timedelta + d0 = datetime(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) + d1 = datetime(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) + date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() + date_dt = np.float64(date_dt_base) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") - else: - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%d") - master_dt = master_split[2] + master_time.strftime("T%H:%M:%S") - slave_dt = slave_split[2] + slave_time.strftime("T%H:%M:%S") + date_ct = d0 + (d1 - d0) / 2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + + master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py index e98a6f1..1e7e556 100755 --- a/testautoRIFT_ISCE.py +++ b/testautoRIFT_ISCE.py @@ -815,21 +815,18 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search - from datetime import date, datetime -# d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) -# d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) + from datetime import datetime, timedelta +# d0 = datetime(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) +# d1 = datetime(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) d0 = datetime.strptime(master_dt,"%Y%m%dT%H:%M:%S.%f") d1 = datetime.strptime(slave_dt,"%Y%m%dT%H:%M:%S.%f") - date_dt_base = d1 - d0 - date_dt = np.float64(date_dt_base.days) + date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() + date_dt = np.float64(date_dt_base) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') - else: - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'time_standard_img1':'UTC','absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'flight_direction_img1':flight_direction_m,'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'time_standard_img2':'UTC','absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'flight_direction_img2':flight_direction_s,'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([[0.0356, 0.0501, 0.0266, 0.0622, 0.0357, 0.0501], @@ -871,12 +868,6 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # # master_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+master_MTL_path))[2][1:-2],':') # slave_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+slave_MTL_path))[2][1:-2],':') - master_time = ['00','00','00'] - slave_time = ['00','00','00'] - - from datetime import time as time1 - master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) - slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) import netcdf_output as no pair_type = 'optical' @@ -895,22 +886,19 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import date - d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) - d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(date_dt_base.days) + from datetime import datetime, timedelta + d0 = datetime(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) + d1 = datetime(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) + date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() + date_dt = np.float64(date_dt_base) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") - else: - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%d") - master_dt = master_split[3][0:8] + master_time.strftime("T%H:%M:%S") - slave_dt = slave_split[3][0:8] + slave_time.strftime("T%H:%M:%S") + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + + master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} @@ -947,13 +935,6 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_filename = master_split[0][-3:]+'_'+master_split[2]+'_'+master_split[4][:3]+'_'+os.path.basename(master_path) slave_filename = slave_split[0][-3:]+'_'+slave_split[2]+'_'+slave_split[4][:3]+'_'+os.path.basename(slave_path) - master_time = ['00','00','00'] - slave_time = ['00','00','00'] - - from datetime import time as time1 - master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) - slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) - import netcdf_output as no pair_type = 'optical' detection_method = 'feature' @@ -970,22 +951,19 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - from datetime import date - d0 = date(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) - d1 = date(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(date_dt_base.days) + from datetime import datetime, timedelta + d0 = datetime(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) + d1 = datetime(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) + date_dt_base = (d1 - d0).total_seconds() / timedelta(days=1).total_seconds() + date_dt = np.float64(date_dt_base) if date_dt < 0: raise Exception('Input image 1 must be older than input image 2') - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") - else: - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%d") - master_dt = master_split[2] + master_time.strftime("T%H:%M:%S") - slave_dt = slave_split[2] + slave_time.strftime("T%H:%M:%S") + date_ct = d0 + (d1 - d0) / 2 + date_center = date_ct.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + + master_dt = d0.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') + slave_dt = d1.strftime("%Y%m%dT%H:%M:%S.%f").rstrip('0') IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version}