diff --git a/docs/dir_structure.md b/docs/dir_structure.md index c82ad7a24..c0afb07c1 100644 --- a/docs/dir_structure.md +++ b/docs/dir_structure.md @@ -459,40 +459,43 @@ Below is a recipe to prepare a stack of interferograms from Sentinel-1: + [https://topex.ucsd.edu/gmtsar/tar/sentinel_time_series_2.pdf](https://topex.ucsd.edu/gmtsar/tar/sentinel_time_series_2.pdf) ``` -$DATA_DIR/StHelensEnvDT156 +$DATA_DIR/SanFranBaySenD42 +├── baseline_table.dat +├── supermaster.PRM ├── geometry -│ └── topo_ll.grd +| ├── azimuth_angle.grd +| ├── dem.grd +│ └── incidence_angle.grd ├── interferograms │ ├── 2004114_2004324 -│   │   ├── 20040423.LED -│   │   ├── 20040423.PRM -│   │   ├── 20041119.LED -│   │   ├── 20041119.PRM -│   │   ├── baseline.txt #generated by SAT_baseline -│   │   ├── corr.grd │   │   ├── corr_ll.grd -│   │   ├── unwrap.grd │   │   └── unwrap_ll.grd │   ├── 2004114_2004359 │ └── ... └── mintpy - └── StHelensEnvDT156.txt + └── SanFranBaySenD42.txt ``` The corresponding template options for `load_data`: ```cfg ## manually specify the following attributes since they are missing from gmtsar products +ALOOKS = 8 #[int], number of looks in the azimuth direction +RLOOKS = 32 #[int], number of looks in the range direction HEADING = -168.0 #[float], satellite heading angle, measured from the north in clockwise as positive # One could open the *.kml file in Google Earth and measure it manually ORBIT_DIRECTION = DESCENDING #[ASCENDING, DESCENDING] -mintpy.load.processor = gmtsar +mintpy.load.processor = gmtsar +mintpy.load.metaFile = $DATA_DIR/SanFranBaySenD42/supermaster.PRM +mintpy.load.baselineDir = $DATA_DIR/SanFranBaySenD42/baseline_table.dat ##---------interferogram datasets: -mintpy.load.unwFile = $DATA_DIR/StHelensEnvDT156/interferograms/*/unwrap_ll.grd -mintpy.load.corFile = $DATA_DIR/StHelensEnvDT156/interferograms/*/corr_ll.grd +mintpy.load.unwFile = $DATA_DIR/SanFranBaySenD42/interferograms/*/unwrap_ll*.grd +mintpy.load.corFile = $DATA_DIR/SanFranBaySenD42/interferograms/*/corr_ll*.grd ##---------geometry datasets: -mintpy.load.demFile = $DATA_DIR/StHelensEnvDT156/geometry/topo_ll.grd +mintpy.load.demFile = $DATA_DIR/SanFranBaySenD42/geometry/dem*.grd +mintpy.load.incAngleFile = $DATA_DIR/SanFranBaySenD42/geometry/incidence_angle.grd +mintpy.load.azAngleFile = $DATA_DIR/SanFranBaySenD42/geometry/azimuth_angle.grd ``` ### Gamma ### diff --git a/pyproject.toml b/pyproject.toml index 68bc7bce8..08c06977b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -110,7 +110,7 @@ include-package-data = true dependencies = { file = ["requirements.txt"] } readme = { file = ["docs/README.md"], content-type = "text/markdown" } -# extra requirements: `pip install pysolid[test]` or `pip install .[test]` +# extra requirements: `pip install mintpy[test]` or `pip install .[test]` [tool.setuptools.dynamic.optional-dependencies.test] file = ["tests/requirements.txt"] diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index c49330db9..0452f1637 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -382,7 +382,7 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key): pairsDict = {} for i, dsPath0 in enumerate(dsPathDict[dsName0]): # date string used in the file/dir path - # YYYYDDD for gmtsar [modern Julian date] + # YYYYDDD for gmtsar [day of the year - 1] # YYYYMMDDTHHMM for uavsar # YYYYMMDD for all the others date6s = readfile.read_attribute(dsPath0)['DATE12'].replace('_','-').split('-') diff --git a/src/mintpy/prep_gmtsar.py b/src/mintpy/prep_gmtsar.py index 7a7f54f21..3b8e230fa 100644 --- a/src/mintpy/prep_gmtsar.py +++ b/src/mintpy/prep_gmtsar.py @@ -19,6 +19,26 @@ ######################################################################### +def get_ifg_file(ifg_dir, coord_type='geo', fbases=['corr', 'phase', 'phasefilt', 'unwrap']): + ifg_file = None + + # search interferogram data file + for fbase in fbases: + if coord_type == 'geo': + fnames = glob.glob(os.path.join(ifg_dir, f'{fbase}_ll*.grd')) + else: + fnames = glob.glob(os.path.join(ifg_dir, f'{fbase}.grd')) + + if len(fnames) > 0: + ifg_file = fnames[0] + break + + if not ifg_file: + print(f'WARNING: No {coord_type}-coord files found in {ifg_dir} with suffix: {fbases}') + + return ifg_file + + def get_prm_files(ifg_dir): """Get the date1/2.prm files in the interferogram directory.""" prm_files = sorted(glob.glob(os.path.join(ifg_dir, '*.PRM'))) @@ -31,12 +51,10 @@ def get_prm_files(ifg_dir): def get_multilook_number(ifg_dir, fbases=['corr', 'phase', 'phasefilt', 'unwrap']): """Get the number of multilook in range and azimuth direction.""" # grab an arbitrary file in radar-coordiantes - rdr_files = [os.path.join(ifg_dir, f'{i}.grd') for i in fbases] - if len(rdr_files) == 0: - raise ValueError(f'No radar-coord files found in {ifg_dir} with suffix: {fbases}') + rdr_file = get_ifg_file(ifg_dir, coord_type='radar') # read step info into multilook number - ds = gdal.Open(rdr_files[0], gdal.GA_ReadOnly) + ds = gdal.Open(rdr_file, gdal.GA_ReadOnly) transform = ds.GetGeoTransform() rlooks = int(abs(transform[1])) alooks = int(abs(transform[5])) @@ -46,13 +64,12 @@ def get_multilook_number(ifg_dir, fbases=['corr', 'phase', 'phasefilt', 'unwrap' def get_lalo_ref(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasefilt', 'unwrap']): """Get the LAT/LON_REF1/2/3/4 from *_ll.grd file.""" # grab an arbitrary file in geo-coordiantes - geo_files = [os.path.join(ifg_dir, f'{i}_ll.grd') for i in fbases] - if len(geo_files) == 0: - print(f'WARNING: No geo-coord files found in {ifg_dir} with suffix: {fbases}') + geo_file = get_ifg_file(ifg_dir, coord_type='geo') + if not geo_file: return prm_dict # read corners lat/lon info - ds = gdal.Open(geo_files[0], gdal.GA_ReadOnly) + ds = gdal.Open(geo_file, gdal.GA_ReadOnly) transform = ds.GetGeoTransform() x_step = abs(transform[1]) y_step = abs(transform[5]) * -1. @@ -88,13 +105,13 @@ def get_lalo_ref(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasefilt', 'unwra def get_slant_range_distance(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasefilt', 'unwrap']): """Get a constant slant range distance in the image center, for dataset in geo-coord.""" - # grab an arbitrary file in radar-coordiantes - rdr_files = [os.path.join(ifg_dir, f'{i}.grd') for i in fbases] - if len(rdr_files) == 0: + # grab an arbitrary file in geo/radar-coordiantes + geo_file = get_ifg_file(ifg_dir, coord_type='geo') + if not geo_file: raise ValueError(f'No radar-coord files found in {ifg_dir} with suffix: {fbases}') # read width from rdr_file - ds = gdal.Open(rdr_files[0], gdal.GA_ReadOnly) + ds = gdal.Open(geo_file, gdal.GA_ReadOnly) width = ds.RasterXSize near_range = float(prm_dict['STARTING_RANGE']) @@ -105,8 +122,30 @@ def get_slant_range_distance(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasef return prm_dict +def read_baseline_table(fname): + """Read GMTSAR baseline table file. + + Example: baseline_table.dat file + # file_ID yyyyddd.fraction day_cnt b_para b_perp + S1_20141231_ALL_F2 2014364.5885345591 364 -58.814087517826 -16.031404204594 + S1_20150301_ALL_F2 2015059.5885222249 424 -33.273566547687 -17.102628888348 + S1_20150325_ALL_F2 2015083.5885250464 448 -104.664278131966 -142.776284597889 + + Parameters: fname - str, path to the baseline table file + Returns: bdict - dict, perpendicular baseline dictionary + """ + print(f'read baseline time-series from file: {fname}') + fc = np.loadtxt(fname, dtype=str, usecols=(1,4)) + date_list = ptime.yyyyddd2yyyymmdd([x[:7] for x in fc[:,0]]) + bperp_list = fc[:,1].astype(np.float32) + bdict = {} + for date_str, bperp in zip(date_list, bperp_list): + bdict[date_str] = bperp + return bdict + + ######################################################################### -def extract_gmtsar_metadata(unw_file, template_file, rsc_file=None, update_mode=True): +def extract_gmtsar_metadata(meta_file, unw_file, template_file, rsc_file=None, update_mode=True): """Extract metadata from GMTSAR interferogram stack.""" # update_mode: check existing rsc_file @@ -116,22 +155,22 @@ def extract_gmtsar_metadata(unw_file, template_file, rsc_file=None, update_mode= ifg_dir = os.path.dirname(unw_file) # 1. read *.PRM file - prm_file = get_prm_files(ifg_dir)[0] - meta = readfile.read_gmtsar_prm(prm_file) + #prm_file = get_prm_files(ifg_dir)[0] + meta = readfile.read_gmtsar_prm(meta_file) meta['PROCESSOR'] = 'gmtsar' # 2. read template file: HEADING, ORBIT_DIRECTION template = readfile.read_template(template_file) - for key in ['HEADING', 'ORBIT_DIRECTION']: + for key in ['HEADING', 'ORBIT_DIRECTION', 'ALOOKS', 'RLOOKS']: if key in template.keys(): meta[key] = template[key].lower() else: raise ValueError('Attribute {} is missing! Please manually specify it in the template file.') # 3. grab A/RLOOKS from radar-coord data file - meta['ALOOKS'], meta['RLOOKS'] = get_multilook_number(ifg_dir) - meta['AZIMUTH_PIXEL_SIZE'] *= meta['ALOOKS'] - meta['RANGE_PIXEL_SIZE'] *= meta['RLOOKS'] + #meta['ALOOKS'], meta['RLOOKS'] = get_multilook_number(ifg_dir) + meta['AZIMUTH_PIXEL_SIZE'] *= int(meta['ALOOKS']) + meta['RANGE_PIXEL_SIZE'] *= int(meta['RLOOKS']) # 4. grab LAT/LON_REF1/2/3/4 from geo-coord data file meta = get_lalo_ref(ifg_dir, meta) @@ -202,7 +241,7 @@ def prepare_geometry(geom_files, meta, update_mode=True): return -def prepare_stack(unw_files, meta, update_mode=True): +def prepare_stack(unw_files, meta, bperp_dict, update_mode=True): """Prepare .rsc file for all unwrapped interferogram files.""" num_file = len(unw_files) if num_file == 0: @@ -222,22 +261,15 @@ def prepare_stack(unw_files, meta, update_mode=True): ifg_meta.update(readfile.read_gdal_vrt(unw_file)) # add DATE12 - prm_files = get_prm_files(ifg_dir) - date1, date2 = (os.path.splitext(os.path.basename(i))[0] for i in prm_files) + date1, date2 = os.path.basename(ifg_dir).replace('-', '_').split('_') + date1 = ptime.yyyyddd2yyyymmdd(date1) + date2 = ptime.yyyyddd2yyyymmdd(date2) ifg_meta['DATE12'] = f'{ptime.yymmdd(date1)}-{ptime.yymmdd(date2)}' # and P_BASELINE_TOP/BOTTOM_HDR - baseline_file = os.path.join(ifg_dir, 'baseline.txt') - if os.path.isfile(baseline_file): - bDict = readfile.read_template(baseline_file) - ifg_meta['P_BASELINE_TOP_HDR'] = bDict['B_perpendicular'] - ifg_meta['P_BASELINE_BOTTOM_HDR'] = bDict['B_perpendicular'] - else: - ifg_meta['P_BASELINE_TOP_HDR'] = '0' - ifg_meta['P_BASELINE_BOTTOM_HDR'] = '0' - msg = f'WARNING: No baseline file found in: {baseline_file}. ' - msg += 'Set P_BASELINE* to 0 and continue.' - print(msg) + bperp = bperp_dict[date2] - bperp_dict[date1] + ifg_meta['P_BASELINE_TOP_HDR'] = bperp + ifg_meta['P_BASELINE_BOTTOM_HDR'] = bperp # write .rsc file rsc_file = unw_file+'.rsc' @@ -256,22 +288,30 @@ def prepare_stack(unw_files, meta, update_mode=True): def prep_gmtsar(inps): # read file path from template file template = readfile.read_template(inps.template_file) + inps.meta_file = glob.glob(template['mintpy.load.metaFile'])[0] + inps.bperp_file = glob.glob(template['mintpy.load.baselineDir'])[0] inps.unw_files = sorted(glob.glob(template['mintpy.load.unwFile'])) inps.cor_files = sorted(glob.glob(template['mintpy.load.corFile'])) inps.dem_file = glob.glob(template['mintpy.load.demFile'])[0] + inps.inc_angle_file = glob.glob(template['mintpy.load.incAngleFile'])[0] + inps.az_angle_file = glob.glob(template['mintpy.load.azAngleFile'])[0] # extract common metadata - rsc_file = os.path.join(inps.mintpy_dir, 'inputs/data.rsc') + rsc_file = os.path.join(os.path.dirname(inps.meta_file), 'data.rsc') meta = extract_gmtsar_metadata( + meta_file=inps.meta_file, unw_file=inps.unw_files[0], template_file=inps.template_file, rsc_file=rsc_file, update_mode=inps.update_mode, ) + # read baseline + bperp_dict = read_baseline_table(inps.bperp_file) + # prepare metadata for geometry files prepare_geometry( - geom_files=[inps.dem_file], + geom_files=[inps.dem_file, inps.inc_angle_file, inps.az_angle_file], meta=meta, update_mode=inps.update_mode, ) @@ -280,6 +320,7 @@ def prep_gmtsar(inps): prepare_stack( unw_files=inps.unw_files, meta=meta, + bperp_dict=bperp_dict, update_mode=inps.update_mode, ) diff --git a/src/mintpy/utils/ptime.py b/src/mintpy/utils/ptime.py index 09b3196ee..d0b8cb631 100644 --- a/src/mintpy/utils/ptime.py +++ b/src/mintpy/utils/ptime.py @@ -236,6 +236,29 @@ def yymmdd2yyyymmdd(date): return date +def yyyyddd2yyyymmdd(date_in): + """Convert GMTSAR folder name in YYYYDDD into YYYYMMDD. + Parameters: date_in - str/list, GMTSAR date format in YYYYDDD, where DDD is the day of year - 1 + Returns: date_out - str/list, date in YYYYMMDD format + """ + if isinstance(date_in, str): + year, doy = date_in[:4], date_in[-3:] + dt_obj = dt.datetime(int(year), 1, 1) + dt.timedelta(days=int(doy)) + date_out = dt_obj.strftime('%Y%m%d') + + elif isinstance(date_in, list): + date_out = [] + for date_str in date_in: + year, doy = date_str[:4], date_str[-3:] + dt_obj = dt.datetime(int(year), 1, 1) + dt.timedelta(days=int(doy)) + date_out.append(dt_obj.strftime('%Y%m%d')) + + else: + return None + + return date_out + + def yy2yyyy(year): """Convert year str from YY to YYYY format""" if year[0] == '9': diff --git a/src/mintpy/utils/readfile.py b/src/mintpy/utils/readfile.py index 6f0f1cd89..39a406310 100644 --- a/src/mintpy/utils/readfile.py +++ b/src/mintpy/utils/readfile.py @@ -207,6 +207,21 @@ 'nan' : np.nan, } +GMTSAR_SENSOR_ID2NAME = { + 1 : 'ers', + 2 : 'ers', + 3 : 'rs1', + 4 : 'env', + 5 : 'alos', + 6 : 'tsx', + 7 : 'csk', + 8 : 'csk', + 9 : 'rs2', + 10: 'sen', + 11: 'gf3', + 12: 'lt1', +} + ######################################################################### def numpy_to_gdal_dtype(np_dtype: DTypeLike) -> int: @@ -1816,6 +1831,9 @@ def _attribute_gmtsar2roipac(prm_dict_in): t_center = dt_center - int(dt_center) prm_dict['CENTER_LINE_UTC'] = str(t_center * 24. * 60. * 60.) + # SC_identity -> PLATFORM + prm_dict['PLATFORM'] = GMTSAR_SENSOR_ID2NAME[int(prm_dict['SC_identity'])] + return prm_dict