From 02b6701aa2c4bc1848277413a7feae8656ff9c15 Mon Sep 17 00:00:00 2001 From: Yang Lei Date: Thu, 28 Oct 2021 09:50:39 -0700 Subject: [PATCH] update support for S2 filenaming conventions and remove projected velocity products --- netcdf_output.py | 400 +++++++++++++++++++++---------------------- testautoRIFT.py | 17 +- testautoRIFT_ISCE.py | 17 +- 3 files changed, 228 insertions(+), 206 deletions(-) diff --git a/netcdf_output.py b/netcdf_output.py index e398c1a..e9859f3 100755 --- a/netcdf_output.py +++ b/netcdf_output.py @@ -1080,206 +1080,206 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SSM1 - varname='vxp' - datatype=np.dtype('int16') - dimensions=('y','x') - FillValue=NoDataValue - var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) - - - var.setncattr('standard_name','projected_x_velocity') - var.setncattr('description','x-direction velocity determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected vx estimates are used') - var.setncattr('units','m/y') - - - if stable_count_p != 0: - 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() - 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))]) - else: - vxp_error_slow = np.nan - if stable_shift_applied_p == 1: - vxp_error = vxp_error_mask - elif stable_shift_applied_p == 2: - vxp_error = vxp_error_slow - else: - vxp_error = vxp_error_mod - var.setncattr('vxp_error',int(round(vxp_error*10))/10) - var.setncattr('vxp_error_description','best estimate of projected_x_velocity error: vxp_error is populated according to the approach used for the velocity bias correction as indicated in "flag_stable_shift"') - - - if stable_shift_applied_p == 2: - var.setncattr('stable_shift',int(round(vxp_mean_shift1*10))/10) - elif stable_shift_applied_p == 1: - var.setncattr('stable_shift',int(round(vxp_mean_shift*10))/10) - else: - var.setncattr('stable_shift',np.nan) - var.setncattr('flag_stable_shift',stable_shift_applied_p) - var.setncattr('flag_stable_shift_description','flag for applying velocity bias correction: 0 = no correction; 1 = correction from overlapping stable surface mask (stationary or slow-flowing surfaces with velocity < 15 m/yr)(top priority); 2 = correction from slowest 25% of overlapping velocities (second priority)') - - - var.setncattr('stable_count_mask',stable_count_p) - var.setncattr('stable_count_slow',stable_count1_p) - if stable_count_p != 0: - var.setncattr('stable_shift_mask',int(round(vxp_mean_shift*10))/10) - else: - var.setncattr('stable_shift_mask',np.nan) - if stable_count1_p != 0: - var.setncattr('stable_shift_slow',int(round(vxp_mean_shift1*10))/10) - else: - var.setncattr('stable_shift_slow',np.nan) - - if stable_count_p != 0: - var.setncattr('vxp_error_mask',int(round(vxp_error_mask*10))/10) - else: - var.setncattr('vxp_error_mask',np.nan) - var.setncattr('vxp_error_mask_description','RMSE over stable surfaces, stationary or slow-flowing surfaces with velocity < 15 m/yr identified from an external mask') - if stable_count1_p != 0: - var.setncattr('vxp_error_slow',int(round(vxp_error_slow*10))/10) - else: - var.setncattr('vxp_error_slow',np.nan) - var.setncattr('vxp_error_slow_description','RMSE over slowest 25% of retrieved velocities') - var.setncattr('vxp_error_modeled',int(round(vxp_error_mod*10))/10) - var.setncattr('vxp_error_modeled_description','1-sigma error calculated using a modeled error-dt relationship') - - - var.setncattr('grid_mapping',mapping_name) - - VXP[noDataMask] = NoDataValue - var[:] = np.round(np.clip(VXP, -32768, 32767)).astype(np.int16) -# var.setncattr('missing_value',np.int16(NoDataValue)) - - - - - - - varname='vyp' - datatype=np.dtype('int16') - dimensions=('y','x') - FillValue=NoDataValue - var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) - - - var.setncattr('standard_name','projected_y_velocity') - var.setncattr('description','y-direction velocity determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected vy estimates are used') - var.setncattr('units','m/y') - - - if stable_count_p != 0: - 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() - 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))]) - else: - vyp_error_slow = np.nan - if stable_shift_applied_p == 1: - vyp_error = vyp_error_mask - elif stable_shift_applied_p == 2: - vyp_error = vyp_error_slow - else: - vyp_error = vyp_error_mod - var.setncattr('vyp_error',int(round(vyp_error*10))/10) - var.setncattr('vyp_error_description','best estimate of projected_y_velocity error: vyp_error is populated according to the approach used for the velocity bias correction as indicated in "flag_stable_shift"') - - - if stable_shift_applied_p == 2: - var.setncattr('stable_shift',int(round(vyp_mean_shift1*10))/10) - elif stable_shift_applied_p == 1: - var.setncattr('stable_shift',int(round(vyp_mean_shift*10))/10) - else: - var.setncattr('stable_shift',np.nan) - var.setncattr('flag_stable_shift',stable_shift_applied_p) - var.setncattr('flag_stable_shift_description','flag for applying velocity bias correction: 0 = no correction; 1 = correction from overlapping stable surface mask (stationary or slow-flowing surfaces with velocity < 15 m/yr)(top priority); 2 = correction from slowest 25% of overlapping velocities (second priority)') - - - var.setncattr('stable_count_mask',stable_count_p) - var.setncattr('stable_count_slow',stable_count1_p) - if stable_count_p != 0: - var.setncattr('stable_shift_mask',int(round(vyp_mean_shift*10))/10) - else: - var.setncattr('stable_shift_mask',np.nan) - if stable_count1_p != 0: - var.setncattr('stable_shift_slow',int(round(vyp_mean_shift1*10))/10) - else: - var.setncattr('stable_shift_slow',np.nan) - - if stable_count_p != 0: - var.setncattr('vyp_error_mask',int(round(vyp_error_mask*10))/10) - else: - var.setncattr('vyp_error_mask',np.nan) - var.setncattr('vyp_error_mask_description','RMSE over stable surfaces, stationary or slow-flowing surfaces with velocity < 15 m/yr identified from an external mask') - if stable_count1_p != 0: - var.setncattr('vyp_error_slow',int(round(vyp_error_slow*10))/10) - else: - var.setncattr('vyp_error_slow',np.nan) - var.setncattr('vyp_error_slow_description','RMSE over slowest 25% of retrieved velocities') - var.setncattr('vyp_error_modeled',int(round(vyp_error_mod*10))/10) - var.setncattr('vyp_error_modeled_description','1-sigma error calculated using a modeled error-dt relationship') - - - var.setncattr('grid_mapping',mapping_name) - - VYP[noDataMask] = NoDataValue - var[:] = np.round(np.clip(VYP, -32768, 32767)).astype(np.int16) -# var.setncattr('missing_value',np.int16(NoDataValue)) - - - - - - - varname='vp' - datatype=np.dtype('int16') - dimensions=('y','x') - FillValue=NoDataValue - var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) - var.setncattr('standard_name','projected_velocity') - var.setncattr('description','velocity magnitude determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected v estimates are used') - var.setncattr('units','m/y') - - var.setncattr('grid_mapping',mapping_name) - - VP[noDataMask] = NoDataValue - var[:] = np.round(np.clip(VP, -32768, 32767)).astype(np.int16) -# var.setncattr('missing_value',np.int16(NoDataValue)) - - - - - - vp_error = v_error_cal(vxp_error, vyp_error) - varname='vp_error' - datatype=np.dtype('int16') - dimensions=('y','x') - FillValue=NoDataValue - var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) - var.setncattr('standard_name','projected_velocity_error') - var.setncattr('description','velocity magnitude error determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected v_error estimates are used') - var.setncattr('units','m/y') - - var.setncattr('grid_mapping',mapping_name) - - VP_error = np.sqrt((vxp_error * VXP / VP)**2 + (vyp_error * VYP / VP)**2) - VP_error[VP==0] = vp_error - VP_error[noDataMask] = NoDataValue - var[:] = np.round(np.clip(VP_error, -32768, 32767)).astype(np.int16) -# var.setncattr('missing_value',np.int16(NoDataValue)) +# varname='vxp' +# datatype=np.dtype('int16') +# dimensions=('y','x') +# FillValue=NoDataValue +# var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) +# +# +# var.setncattr('standard_name','projected_x_velocity') +# var.setncattr('description','x-direction velocity determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected vx estimates are used') +# var.setncattr('units','m/y') +# +# +# if stable_count_p != 0: +# 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() - 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))]) +# else: +# vxp_error_slow = np.nan +# if stable_shift_applied_p == 1: +# vxp_error = vxp_error_mask +# elif stable_shift_applied_p == 2: +# vxp_error = vxp_error_slow +# else: +# vxp_error = vxp_error_mod +# var.setncattr('vxp_error',int(round(vxp_error*10))/10) +# var.setncattr('vxp_error_description','best estimate of projected_x_velocity error: vxp_error is populated according to the approach used for the velocity bias correction as indicated in "flag_stable_shift"') +# +# +# if stable_shift_applied_p == 2: +# var.setncattr('stable_shift',int(round(vxp_mean_shift1*10))/10) +# elif stable_shift_applied_p == 1: +# var.setncattr('stable_shift',int(round(vxp_mean_shift*10))/10) +# else: +# var.setncattr('stable_shift',np.nan) +# var.setncattr('flag_stable_shift',stable_shift_applied_p) +# var.setncattr('flag_stable_shift_description','flag for applying velocity bias correction: 0 = no correction; 1 = correction from overlapping stable surface mask (stationary or slow-flowing surfaces with velocity < 15 m/yr)(top priority); 2 = correction from slowest 25% of overlapping velocities (second priority)') +# +# +# var.setncattr('stable_count_mask',stable_count_p) +# var.setncattr('stable_count_slow',stable_count1_p) +# if stable_count_p != 0: +# var.setncattr('stable_shift_mask',int(round(vxp_mean_shift*10))/10) +# else: +# var.setncattr('stable_shift_mask',np.nan) +# if stable_count1_p != 0: +# var.setncattr('stable_shift_slow',int(round(vxp_mean_shift1*10))/10) +# else: +# var.setncattr('stable_shift_slow',np.nan) +# +# if stable_count_p != 0: +# var.setncattr('vxp_error_mask',int(round(vxp_error_mask*10))/10) +# else: +# var.setncattr('vxp_error_mask',np.nan) +# var.setncattr('vxp_error_mask_description','RMSE over stable surfaces, stationary or slow-flowing surfaces with velocity < 15 m/yr identified from an external mask') +# if stable_count1_p != 0: +# var.setncattr('vxp_error_slow',int(round(vxp_error_slow*10))/10) +# else: +# var.setncattr('vxp_error_slow',np.nan) +# var.setncattr('vxp_error_slow_description','RMSE over slowest 25% of retrieved velocities') +# var.setncattr('vxp_error_modeled',int(round(vxp_error_mod*10))/10) +# var.setncattr('vxp_error_modeled_description','1-sigma error calculated using a modeled error-dt relationship') +# +# +# var.setncattr('grid_mapping',mapping_name) +# +# VXP[noDataMask] = NoDataValue +# var[:] = np.round(np.clip(VXP, -32768, 32767)).astype(np.int16) +## var.setncattr('missing_value',np.int16(NoDataValue)) +# +# +# +# +# +# +# varname='vyp' +# datatype=np.dtype('int16') +# dimensions=('y','x') +# FillValue=NoDataValue +# var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) +# +# +# var.setncattr('standard_name','projected_y_velocity') +# var.setncattr('description','y-direction velocity determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected vy estimates are used') +# var.setncattr('units','m/y') +# +# +# if stable_count_p != 0: +# 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() - 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))]) +# else: +# vyp_error_slow = np.nan +# if stable_shift_applied_p == 1: +# vyp_error = vyp_error_mask +# elif stable_shift_applied_p == 2: +# vyp_error = vyp_error_slow +# else: +# vyp_error = vyp_error_mod +# var.setncattr('vyp_error',int(round(vyp_error*10))/10) +# var.setncattr('vyp_error_description','best estimate of projected_y_velocity error: vyp_error is populated according to the approach used for the velocity bias correction as indicated in "flag_stable_shift"') +# +# +# if stable_shift_applied_p == 2: +# var.setncattr('stable_shift',int(round(vyp_mean_shift1*10))/10) +# elif stable_shift_applied_p == 1: +# var.setncattr('stable_shift',int(round(vyp_mean_shift*10))/10) +# else: +# var.setncattr('stable_shift',np.nan) +# var.setncattr('flag_stable_shift',stable_shift_applied_p) +# var.setncattr('flag_stable_shift_description','flag for applying velocity bias correction: 0 = no correction; 1 = correction from overlapping stable surface mask (stationary or slow-flowing surfaces with velocity < 15 m/yr)(top priority); 2 = correction from slowest 25% of overlapping velocities (second priority)') +# +# +# var.setncattr('stable_count_mask',stable_count_p) +# var.setncattr('stable_count_slow',stable_count1_p) +# if stable_count_p != 0: +# var.setncattr('stable_shift_mask',int(round(vyp_mean_shift*10))/10) +# else: +# var.setncattr('stable_shift_mask',np.nan) +# if stable_count1_p != 0: +# var.setncattr('stable_shift_slow',int(round(vyp_mean_shift1*10))/10) +# else: +# var.setncattr('stable_shift_slow',np.nan) +# +# if stable_count_p != 0: +# var.setncattr('vyp_error_mask',int(round(vyp_error_mask*10))/10) +# else: +# var.setncattr('vyp_error_mask',np.nan) +# var.setncattr('vyp_error_mask_description','RMSE over stable surfaces, stationary or slow-flowing surfaces with velocity < 15 m/yr identified from an external mask') +# if stable_count1_p != 0: +# var.setncattr('vyp_error_slow',int(round(vyp_error_slow*10))/10) +# else: +# var.setncattr('vyp_error_slow',np.nan) +# var.setncattr('vyp_error_slow_description','RMSE over slowest 25% of retrieved velocities') +# var.setncattr('vyp_error_modeled',int(round(vyp_error_mod*10))/10) +# var.setncattr('vyp_error_modeled_description','1-sigma error calculated using a modeled error-dt relationship') +# +# +# var.setncattr('grid_mapping',mapping_name) +# +# VYP[noDataMask] = NoDataValue +# var[:] = np.round(np.clip(VYP, -32768, 32767)).astype(np.int16) +## var.setncattr('missing_value',np.int16(NoDataValue)) +# +# +# +# +# +# +# varname='vp' +# datatype=np.dtype('int16') +# dimensions=('y','x') +# FillValue=NoDataValue +# var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) +# var.setncattr('standard_name','projected_velocity') +# var.setncattr('description','velocity magnitude determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected v estimates are used') +# var.setncattr('units','m/y') +# +# var.setncattr('grid_mapping',mapping_name) +# +# VP[noDataMask] = NoDataValue +# var[:] = np.round(np.clip(VP, -32768, 32767)).astype(np.int16) +## var.setncattr('missing_value',np.int16(NoDataValue)) +# +# +# +# +# +# vp_error = v_error_cal(vxp_error, vyp_error) +# varname='vp_error' +# datatype=np.dtype('int16') +# dimensions=('y','x') +# FillValue=NoDataValue +# var = nc_outfile.createVariable(varname,datatype,dimensions, fill_value=FillValue, zlib=True, complevel=2, shuffle=True, chunksizes=ChunkSize) +# var.setncattr('standard_name','projected_velocity_error') +# var.setncattr('description','velocity magnitude error determined by projecting radar range measurements onto an a priori flow vector. Where projected errors are larger than those determined from range and azimuth measurements, unprojected v_error estimates are used') +# var.setncattr('units','m/y') +# +# var.setncattr('grid_mapping',mapping_name) +# +# VP_error = np.sqrt((vxp_error * VXP / VP)**2 + (vyp_error * VYP / VP)**2) +# VP_error[VP==0] = vp_error +# VP_error[noDataMask] = NoDataValue +# var[:] = np.round(np.clip(VP_error, -32768, 32767)).astype(np.int16) +## var.setncattr('missing_value',np.int16(NoDataValue)) diff --git a/testautoRIFT.py b/testautoRIFT.py index d799e56..3a5936f 100755 --- a/testautoRIFT.py +++ b/testautoRIFT.py @@ -928,8 +928,19 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_split = master_path.split('_') slave_split = slave_path.split('_') - master_filename = os.path.basename(master_path) - slave_filename = os.path.basename(slave_path) + import re + if re.findall("://",master_path).__len__() > 0: + master_filename_full = master_path.split('/') + for item in master_filename_full: + if re.findall("S2._",item).__len__() > 0: + master_filename = item + slave_filename_full = slave_path.split('/') + for item in slave_filename_full: + if re.findall("S2._",item).__len__() > 0: + slave_filename = item + else: + master_filename = os.path.basename(master_path)[:-8] + slave_filename = os.path.basename(slave_path)[:-8] # 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) @@ -944,7 +955,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search raise Exception('Input search range is all zero everywhere, thus no search conducted') PPP = roi_valid_percentage * 100 if ncname is None: - out_nc_filename = f"./{master_filename[0:-8]}_X_{slave_filename[0:-8]}" \ + out_nc_filename = f"./{master_filename}_X_{slave_filename}" \ f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py index d46f671..25430be 100755 --- a/testautoRIFT_ISCE.py +++ b/testautoRIFT_ISCE.py @@ -932,8 +932,19 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_split = master_path.split('_') slave_split = slave_path.split('_') - master_filename = os.path.basename(master_path) - slave_filename = os.path.basename(slave_path) + import re + if re.findall("://",master_path).__len__() > 0: + master_filename_full = master_path.split('/') + for item in master_filename_full: + if re.findall("S2._",item).__len__() > 0: + master_filename = item + slave_filename_full = slave_path.split('/') + for item in slave_filename_full: + if re.findall("S2._",item).__len__() > 0: + slave_filename = item + else: + master_filename = os.path.basename(master_path)[:-8] + slave_filename = os.path.basename(slave_path)[:-8] # 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) @@ -948,7 +959,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search raise Exception('Input search range is all zero everywhere, thus no search conducted') PPP = roi_valid_percentage * 100 if ncname is None: - out_nc_filename = f"./{master_filename[0:-8]}_X_{slave_filename[0:-8]}" \ + out_nc_filename = f"./{master_filename}_X_{slave_filename}" \ f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc"