@@ -618,6 +618,7 @@ def bistatic_delay(self, range_step=1, az_step=1):
618618
619619 return isce3 .core .LUT2d (slant_vec , az_vec , bistatic_correction )
620620
621+
621622 def geometrical_and_steering_doppler (self , range_step = 500 , az_step = 50 ):
622623 """
623624 Compute total Doppler which is the sum of two components:
@@ -661,6 +662,7 @@ def geometrical_and_steering_doppler(self, range_step=500, az_step=50):
661662
662663 return isce3 .core .LUT2d (range_vec , az_vec , total_doppler )
663664
665+
664666 def doppler_induced_range_shift (self , range_step = 500 , az_step = 50 ):
665667 """
666668 Computes the range delay caused by the Doppler shift as described
@@ -688,6 +690,7 @@ def doppler_induced_range_shift(self, range_step=500, az_step=50):
688690
689691 return isce3 .core .LUT2d (range_vec , az_vec , tau_corr )
690692
693+
691694 def az_carrier_components (self , offset , position ):
692695 '''
693696 Estimate azimuth carrier and store in numpy arrary. Also return
@@ -742,6 +745,7 @@ def az_carrier_components(self, offset, position):
742745
743746 return AzimuthCarrierComponents (kt , eta , eta_ref )
744747
748+
745749 def az_fm_rate_mismatch_mitigation (self ,
746750 path_dem : str ,
747751 path_scratch : str = None ,
@@ -751,10 +755,11 @@ def az_fm_rate_mismatch_mitigation(self,
751755 numiter_rdr2geo = 25 ,
752756 custom_radargrid = None ):
753757 '''
754- Calculate azimuth FM rate mismatch mitigation
755- Based on ETAD-DLR-DD-0008, Algorithm Technical Baseline Document.
756- Available: https://sentinels.copernicus.eu/documents/247904/4629150/
757- Sentinel-1-ETAD-Algorithm-Technical-Baseline-Document.pdf
758+ - Calculate Lon / Lat / Hgt in radar grid, to be used for the
759+ actual computation of az fm mismatch rate
760+ - Define the radar grid for the correction.
761+ - call `az_fm_rate_mismatch_from_llh` to do the actual computation
762+
758763
759764 Parameters:
760765 -----------
@@ -779,9 +784,10 @@ def az_fm_rate_mismatch_mitigation(self,
779784 -------
780785 _: isce3.core.LUT2d
781786 azimuth FM rate mismatch rate in radar grid in seconds.
782- Examples
783- ----------
784- >>> correction_grid = burst.az_fm_rate_mismatch_mitigation("my_dem.tif")
787+
788+ Examples
789+ ----------
790+ >>> correction_grid = burst.az_fm_rate_mismatch_mitigation("my_dem.tif")
785791 '''
786792
787793 # Create temporary directory for scratch when
@@ -792,9 +798,9 @@ def az_fm_rate_mismatch_mitigation(self,
792798 else :
793799 temp_dir_obj = None
794800
795-
796801 os .makedirs (path_scratch , exist_ok = True )
797802
803+ # define the radar grid to calculate az fm mismatch rate
798804 correction_radargrid = self .as_isce3_radargrid ()
799805 if custom_radargrid is not None :
800806 correction_radargrid = custom_radargrid
@@ -819,6 +825,97 @@ def az_fm_rate_mismatch_mitigation(self,
819825 width_radargrid ,
820826 self .as_isce3_radargrid ().ref_epoch )
821827
828+ # Run topo on scratch directory
829+ if not os .path .isfile (path_dem ):
830+ raise FileNotFoundError (f'not found - DEM { path_dem } )
831+ dem_raster = isce3 .io .Raster (path_dem )
832+ epsg = dem_raster .get_epsg ()
833+ proj = isce3 .core .make_projection (epsg )
834+ ellipsoid = proj .ellipsoid
835+ grid_doppler = isce3 .core .LUT2d ()
836+
837+ rdr2geo_obj = isce3 .geometry .Rdr2Geo (
838+ correction_radargrid ,
839+ self .orbit ,
840+ ellipsoid ,
841+ grid_doppler ,
842+ threshold = threshold_rdr2geo ,
843+ numiter = numiter_rdr2geo )
844+
845+ str_datetime = datetime .datetime .now ().strftime ('%Y%m%d_%H%M%S%f' )
846+ list_filename_llh = [f'{ path_scratch } { llh } { str_datetime }
847+ for llh
848+ in ['lat' , 'lon' , 'hgt' ]]
849+ lat_raster , lon_raster , hgt_raster = [isce3 .io .Raster (
850+ filename_llh ,
851+ correction_radargrid .width ,
852+ correction_radargrid .length ,
853+ 1 ,
854+ gdal .GDT_Float64 ,
855+ 'ENVI' )
856+ for filename_llh in list_filename_llh ]
857+
858+ rdr2geo_obj .topo (dem_raster , lon_raster , lat_raster , hgt_raster )
859+
860+ # make sure that the ISCE3 rasters are written out to file system
861+ lat_raster .close_dataset ()
862+ lon_raster .close_dataset ()
863+ hgt_raster .close_dataset ()
864+
865+ # Load the lon / lat / hgt value from the raster
866+ arr_lat = gdal .Open (list_filename_llh [0 ], gdal .GA_ReadOnly ).ReadAsArray ()
867+ arr_lon = gdal .Open (list_filename_llh [1 ], gdal .GA_ReadOnly ).ReadAsArray ()
868+ arr_hgt = gdal .Open (list_filename_llh [2 ], gdal .GA_ReadOnly ).ReadAsArray ()
869+
870+ lut_az_fm_mismatch = self .az_fm_rate_mismatch_from_llh (arr_lat ,
871+ arr_lon ,
872+ arr_hgt ,
873+ ellipsoid ,
874+ correction_radargrid )
875+
876+ # Clean up the temporary directory in case it exists.
877+ if temp_dir_obj is not None :
878+ temp_dir_obj .cleanup ()
879+
880+ return lut_az_fm_mismatch
881+
882+
883+ def az_fm_rate_mismatch_from_llh (self ,
884+ lat_map : np .ndarray ,
885+ lon_map : np .ndarray ,
886+ hgt_map : np .ndarray ,
887+ ellipsoid : isce3 .core .Ellipsoid ,
888+ correction_radargrid : isce3 .product .RadarGridParameters ):
889+ '''
890+ Take in lat / lon / hgt in radar grid, along with ellipsoid and the radar grid.
891+ Calculate azimuth FM rate mismatch mitigation based on algorithm
892+ described in [1]
893+
894+ Parameters:
895+ -----------
896+ lat_map: np.ndarray
897+ Latitude in radar grid, unit: degrees
898+ lon_map: np.ndarray,
899+ Longitude in radar grid, unit: degrees
900+ hgt_map: np.ndarray,
901+ Height in radar grid, unit: meters
902+ ellipsoid: isce3.core.Ellipsoid
903+ Reference ellipsoid to be used
904+ correction_radargrid: isce3.product.RadarGridParameters
905+ Radar grid as the definition of the correction grid
906+
907+ Return:
908+ -------
909+ _: isce3.core.LUT2d
910+ azimuth FM rate mismatch rate in radar grid in seconds.
911+
912+ References
913+ ----------
914+ [1] C. Gisinger, "S-1 ETAD project Algorithm Technical
915+ Baseline Document, The integration of GIS, remote sensing,"
916+ DLR, ETAD-DLR-DD-0008, 2020.
917+ '''
918+
822919 # Define the correction grid from the radargrid
823920 # Also define the staggered grid in azimuth to calculate acceeration
824921 width_grid = correction_radargrid .width
@@ -842,6 +939,7 @@ def az_fm_rate_mismatch_mitigation(self,
842939 vec_position_intp = np .zeros ((length_grid , 3 ))
843940 vec_vel_intp = np .zeros ((length_grid , 3 ))
844941 vec_vel_intp_staggered = [None ] * (length_grid + 1 )
942+
845943 for i_azimuth , t_azimuth in enumerate (vec_t ):
846944 vec_position_intp [i_azimuth , :], vec_vel_intp [i_azimuth , :] = \
847945 self .orbit .interpolate (t_azimuth )
@@ -880,7 +978,7 @@ def az_fm_rate_mismatch_mitigation(self,
880978
881979 # Interpolate the DC and fm rate coeffs along azimuth time
882980 def interp_coeffs (az_time , coeffs , az_time_interp ):
883- '''Convenience function to interpolate DC and FM rate coeffiicients
981+ '''Convenience function to interpolate DC and FM rate coefficients
884982 '''
885983 interpolated_coeffs = []
886984 for i in range (3 ):
@@ -897,41 +995,6 @@ def interp_coeffs(az_time, coeffs, az_time_interp):
897995 self .extended_coeffs .fm_rate_coeff_arr ,
898996 vec_t )
899997
900- # Run topo on scratch directory
901- dem_raster = isce3 .io .Raster (path_dem )
902- epsg = dem_raster .get_epsg ()
903- proj = isce3 .core .make_projection (epsg )
904- ellipsoid = proj .ellipsoid
905- grid_doppler = isce3 .core .LUT2d ()
906-
907- rdr2geo_obj = isce3 .geometry .Rdr2Geo (
908- correction_radargrid ,
909- self .orbit ,
910- ellipsoid ,
911- grid_doppler ,
912- threshold = threshold_rdr2geo ,
913- numiter = numiter_rdr2geo )
914-
915- str_datetime = datetime .datetime .now ().strftime ('%Y%m%d_%H%M%S%f' )
916- list_filename_llh = [f'{ path_scratch } { llh } { str_datetime }
917- for llh
918- in ['lat' , 'lon' , 'hgt' ]]
919- lat_raster , lon_raster , hgt_raster = [isce3 .io .Raster (
920- filename_llh ,
921- correction_radargrid .width ,
922- correction_radargrid .length ,
923- 1 ,
924- gdal .GDT_Float64 ,
925- 'ENVI' )
926- for filename_llh in list_filename_llh ]
927-
928- rdr2geo_obj .topo (dem_raster , lon_raster , lat_raster , hgt_raster )
929-
930- # make sure that the ISCE3 rasters are written out to file system
931- lat_raster .close_dataset ()
932- lon_raster .close_dataset ()
933- hgt_raster .close_dataset ()
934-
935998 # Do the bunch of computation
936999 kappa_steer_vec = (2 * (np .linalg .norm (vec_vel_intp , axis = 1 ))
9371000 / isce3 .core .speed_of_light * self .radar_center_frequency
@@ -973,18 +1036,10 @@ def interp_coeffs(az_time, coeffs, az_time_interp):
9731036 + (freq_dcg_burst / kappa_annotation_burst
9741037 - grid_freq_dcg / kappa_annotation_grid )))
9751038
976- lat_map = gdal .Open (list_filename_llh [0 ], gdal .GA_ReadOnly ).ReadAsArray ()
977- lon_map = gdal .Open (list_filename_llh [1 ], gdal .GA_ReadOnly ).ReadAsArray ()
978- hgt_map = gdal .Open (list_filename_llh [2 ], gdal .GA_ReadOnly ).ReadAsArray ()
979-
980- # Clean up the temporary directory in case it exists.
981- if temp_dir_obj is not None :
982- temp_dir_obj .cleanup ()
983-
9841039 x_ecef , y_ecef , z_ecef = _llh_to_ecef (lat_map ,
985- lon_map ,
986- hgt_map ,
987- ellipsoid )
1040+ lon_map ,
1041+ hgt_map ,
1042+ ellipsoid )
9881043
9891044 # Populate the position, velocity, and
9901045 # acceleration of the satellite to the correction grid
0 commit comments