===================================================== Cassini RADAR Titan Plains Backscatter Analysis ===================================================== OVERVIEW -------- This repository contains the MATLAB code and processed data used for the analysis presented in the paper "Titan's Plains Revealed: Evidence for a Layered Surface" (Fine et al., JGR: Planets). The workflow processes Cassini RADAR data (Synthetic Aperture Radar - SAR, and Altimetry) for selected undifferentiated plains regions on Titan, coregisters multiple SAR observations, extracts backscatter (sigma0) and incidence angle data, fits various backscatter models (including single-surface and two-layer models), and generates figures presented in the paper. The core finding is that a two-layer model, consisting of a thin (approx < 1m), low-density, smooth organic layer overlying a denser, rougher substrate, is required to accurately describe the multi-angle radar backscatter from Titan's plains. PREREQUISITES ------------- - MATLAB: The code was developed and tested using MATLAB (version R2023b recommended). - MATLAB Toolboxes: - Image Processing Toolbox (required for Coreg.m and image manipulation) - Statistics and Machine Learning Toolbox (required for fitting, prctile, etc.) - Mapping Toolbox (may be implicitly used by some functions or for visualization) - Cassini RADAR Data: Raw or calibrated Cassini RADAR data (BIDR for SAR, LBDR/ABDR for Altimetry) are needed as input if regenerating site data from scratch. These are available from the NASA Planetary Data System (PDS) Imaging Node: https://pds-imaging.jpl.nasa.gov/volumes/radar.html - Specific BIS*.mat files (calibrated, noise-subtracted SAR data) are used by SiteAnalysis.m. - Specific ABDR_SUMMARY*.csv files (Altimeter Burst Data Record summaries) are used by SiteAnalysis.m. - Specific TXXX_burstID_sigma0.mat files (reprocessed altimetry backscatter) are used by SiteAnalysis.m. DATA DESCRIPTION ---------------- - Raw Input Data (Not included, see Prerequisites): Calibrated SAR images (.mat format derived from PDS BIDRs) and Altimetry summary files (.csv format derived from PDS ABDRs) or reprocessed altimetry .mat files. - Processed Site Data (Site[A-M].mat): These files contain the processed data for each individual study site, generated by SiteAnalysis.m. Each .mat file typically includes: - alt_inc, alt_sigma0: Altimetry incidence angles (degrees) and backscatter (linear power units) for the site, if available. - inc, sigma0: SAR incidence angles (degrees) and backscatter (linear power units) for each coregistered SAR track, extracted from the masked ROI. Dimensions are (N pixels in ROI, M tracks). - co_inc, co_sar: Coregistered SAR incidence angle and backscatter image cubes (pixels x pixels x M tracks) before masking. - mask_inc, mask_sar: Coregistered SAR image cubes, masked to the selected ROI (NaN outside ROI). - mask: Logical mask defining the ROI within the coregistered image dimensions. - sar_filenames: List of input SAR filenames for the site. - site: Name of the site (e.g., 'SiteA'). - data: Structure containing combined SAR and Altimetry data (sigma0 [linear], inc [radians]) for the site, split into data.alt (with altimetry) and data.noalt (SAR only). Note: This structure is primarily for intermediate use within SiteAnalysis.m. - Repository Data: The Site[A-M].mat files used for the final analysis in the paper are included in this repository. WORKFLOW OVERVIEW ----------------- The analysis proceeds in two main stages: individual site processing and aggregate analysis/plotting. Individual Site Processing (SiteAnalysis.m): ------------------------------------------ - Define Site: Select a site section within the script. Define input SAR filenames (sar_filenames) and altimetry information (csv filename and relevant burst IDs ind). - Coregistration (Coreg.m): Run the Coreg section. Interactively crop the primary SAR image, select tie points on features in subsequent images, and refine tie points for precise coregistration of all SAR tracks for the site. - Masking (Mask.m): Run the Mask section. Interactively draw a polygon on the coregistered images to define the Region of Interest (ROI), ensuring homogeneity and avoiding edge artifacts across all tracks. - Data Extraction: Run the subsequent sections to extract sigma0 and inc data from the masked ROI for each track, combine with altimetry data (if applicable), and format it. - (Optional) Individual Fit: Run the fitBackscatter section to perform model fits on the data from this single site. - Save Data: Manually save the relevant workspace variables (alt_inc, alt_sigma0, inc, sigma0, co_inc, co_sar, mask_inc, mask_sar, mask, sar_filenames, site) into a SiteX.mat file (e.g., SiteA.mat). --- Aggregate Analysis and Plotting (allSites.m): ------------------------------------------- - Load Data: Specify the list of SiteX.mat files (site_names) to include in the aggregate analysis. The script loads the inc, sigma0, alt_inc, and alt_sigma0 variables from each specified file. - Combine Data: Concatenate the data from all selected sites into a single dataset (data.alt and data.noalt). - Fit Models (fitBackscatter.m): Fit the defined backscatter models (single-surface and two-layer) to the combined dataset. - Generate Figures: Execute plotting scripts (plotBlackPoints.m, fitsFigure.m, power_phase_plot.m, etc.) to create the figures presented in the paper using the aggregate data and fit results. FILE DESCRIPTIONS ----------------- Main Scripts: - SiteAnalysis.m: Main script for processing data for a single Titan plains site. Includes sections for defining inputs, coregistration, masking, data extraction, and (optionally) fitting. Generates the SiteX.mat files. Requires user interaction. - allSites.m: Main script for the aggregate analysis. Loads multiple SiteX.mat files, combines the data, performs model fitting on the aggregate dataset, and calls plotting functions. Core Processing Functions: - Coreg.m: Function called by SiteAnalysis.m to perform interactive coregistration of multiple SAR images using user-defined tie points. Requires Image Processing Toolbox. - Mask.m: Function called by SiteAnalysis.m for interactive definition of the Region of Interest (ROI) mask on coregistered images. Requires Image Processing Toolbox. - fitBackscatter.m: Performs least-squares fitting of various backscatter models (defined below) to the provided data structure (containing sigma0 and inc). Returns fit objects and parameters. Backscatter Models: - exponential_bscatter.m: Calculates quasi-specular backscatter using the exponential surface model (Eq. 2 in paper). - gaussian_bscatter.m: Calculates quasi-specular backscatter using the Gaussian surface model (Eq. 3 in paper). - hagfors.m: Calculates quasi-specular backscatter using the Hagfors surface model (Eq. 4 in paper). - volume_scatter.m: Calculates the empirical diffuse (volume) scattering component (Eq. 5 in paper). - simple_two_layer.m: Implements the two-layer backscatter model (Eq. 6 in paper), combining surface reflection, subsurface reflection (parameterized by P2), and volume scattering. Note: The parameter fit as eps2 in this function corresponds to the Power Coefficient P2 defined in the paper (Eq. 7), not the subsurface dielectric constant epsilon2. - transmission_function.m: Helper function to calculate the transmission of power from the first layer to the second layer in the two-layer model. Plotting and Figure Generation: - plotFitdB.m: Plots the fitted backscatter curves against the data (binned or raw points) in dB scale. Called by fitBackscatter.m. - plotBlackPoints.m: Generates a plot of all aggregated data points (SAR and Altimetry) used for fitting (similar to points in Fig. 2). Called by allSites.m. - fitsFigure.m: Generates the main model fits figure (Fig. 2 in the paper), showing data and best-fit curves for single-surface and two-layer models. Requires results from fitBackscatter.m run on aggregate data. - power_phase_plot.m: Generates the phase space plot for the two-layer model (Fig. 3 in the paper), showing the relationship between surface thickness, surface loss tangent, and subsurface dielectric constant consistent with the fitted Power Coefficient P2. Note: Requires the P2 value obtained from the simple_two_layer fit (stored as eps2 in the fit object). - contextFigure.m: Generates a draft context map figure (similar Fig. 1 in the paper). - Plotsite.m: Helper function called by SiteAnalysis.m to display the masked, coregistered SAR images for review. - fixfig.m: Utility function for adjusting figure properties for publication quality. - distinguishable_colors.m: Utility function (must be obtained externally) for generating distinct colors for plotting multiple datasets. Utility and Helper Functions: - FindLatLon.m: Utility function to find and display the approximate center latitude/longitude of a site defined in a SiteX.mat file. - findLowInc.m: Utility script to search through SAR metadata to find SAR tracks covering a specific region at desired low incidence angles. Used for initial site identification. - read_isis.m: Utility function for reading specific data formats. - pds_label_parse_v2.m: Utility function for parsing information from PDS label files. - createInfoCsv.m: Utility function to compile metadata info on all of the sites. - calc_chi_square.m: Calculates the chi-square goodness-of-fit statistic for a given model, data, and parameters. - theoretical_Tb.m: Calculates the expected brightness temperature based on the two-layer model parameters using the emissivity model from Le Gall et al. (2016), used for comparison with Cassini radiometry data. Data Files: - SiteA.mat ... SiteM.mat: Processed data files for each study site, generated by SiteAnalysis.m. See Data Description section. USAGE NOTES ----------- - Paths: The scripts may contain hardcoded paths (e.g., in SiteAnalysis.m addpath lines). These will need to be adjusted for your local environment. Running SiteAnalysis.m: ----------------------- 1. Select Site & Inputs: Open SiteAnalysis.m. Choose the desired site section (e.g., %% SITE A DEFINITION). Verify or input the list of SAR BISQ*.mat filenames (sar_filenames). If altimetry is available for the site, provide the corresponding ABDR_SUMMARY*.csv filename (csv) or reprocessed TXXX_burstID_sigma0.mat file and specify the range of relevant Burst IDs (ind). Run this definition section to load the filenames and altimetry info into the workspace. 2. Coregistration (Coreg.m): Run the section starting with [co_sar, co_inc] = Coreg(sar_filenames);. - A window will display the first (primary) SAR track. Interactively draw a rectangle (crop) around the general area of interest and double-click inside it to confirm. - Subsequent SAR tracks for the site will be shown sequentially; crop each one similarly. - After cropping, a zoomed image of the primary track appears. Click on several easily identifiable features (tie points) across the image. Press Enter when done selecting points. - Pairs of images (primary on left, secondary on right) will appear with automatically estimated tie points (red crosses). Click and drag the red crosses on the right image to precisely match the features marked by the blue crosses on the left image. Close the figure window when satisfied with the alignment for that pair. This repeats for all secondary images. The output (co_sar, co_inc) contains the coregistered image stacks. 3. Masking ROI (Mask.m): Run the section starting with [mask_sar, mask_inc, mask] = Mask(co_sar, co_inc);. - Two figures appear. The first is the primary coregistered image. The second (Plotsite) shows all coregistered tracks for comparison. - On the first figure (primary image), draw a polygon outlining the desired final Region of Interest (ROI). Use the zoom tool and refer to the second figure to ensure the chosen area is homogeneous, well-covered by all tracks, and avoids artifacts (like SAR beam edges). Double-click inside the polygon to confirm the mask. - A new figure shows the masked area on the primary image. The outputs mask_sar, mask_inc contain the masked image stacks, and mask is the logical mask itself. 4. Data Extraction & Review: Run the next few sections. These extract the backscatter (sigma0) and incidence angle (inc) values from within the mask for each SAR track, combine them with the specified altimetry data (alt_sigma0, alt_inc), and generate review plots (histograms, scatter plot of backscatter vs. incidence angle). The intermediate structure data (with subfields alt and noalt) is created, containing combined data with angles in radians. 5. (Optional) Fit Single Site: Run the fits = fitBackscatter(data, 1); section to fit the backscatter models to the data from only this site. 6. Save Data: Manually save the essential workspace variables to a .mat file named appropriately (e.g., SiteA.mat). The crucial variables to save are: alt_inc, alt_sigma0, inc, sigma0, co_inc, co_sar, mask_inc, mask_sar, mask, sar_filenames, site. --- Running allSites.m: ------------------- 1. Select Sites: Modify the site_names cell array at the beginning of the script to list the SiteX.mat files you want to include in the combined analysis. 2. Ensure Availability: Make sure the specified SiteX.mat files are in the MATLAB path or the current directory. 3. Run Script: Execute the entire allSites.m script. 4. Process: The script will load inc, sigma0, alt_inc, alt_sigma0 from each specified file, concatenate them into aggregate datasets (data.alt and data.noalt), fit the backscatter models to the aggregate data.alt using fitBackscatter, store results in the fits structure, and generate the main figures (like Fig. 2, Fig. 3) using the aggregate results. 5. Model Parameter Interpretation: Pay close attention to the definition of parameters in the code vs. the paper, especially eps2 in simple_two_layer.m which represents P2 (Power Coefficient), not epsilon2 (subsurface permittivity). The power_phase_plot.m script correctly interprets this fitted value as P2. ALTIMETRY ROUGHNESS ESTIMATION CODE ----------------------------------- This section describes the MATLAB code used for estimating large-scale surface roughness (RMS height, sigma_h) from Cassini RADAR altimetry echo waveforms, specifically for Titan's plains. This corresponds to the analysis described in Section 2.4 and results in Section 3.2 / Table 3 of the paper. Script: - NADIR_ROUGHNESS_EST_PLAINS.m: This script performs the roughness estimation. Functionality: - Loads Data: Prompts the user to select an input `.mat` file (e.g., `T30_WS.mat`). These files are expected to contain variables like: - `REALI_Plains`: Likely the processed altimetry echo waveforms. - `SC_radius_Plains`: Spacecraft altitude/range information. - `teta_off_Plains`: Off-nadir pointing angle for the waveforms. - Preprocessing (Optional): If the `flag` variable is set to 1, the script attempts to realign waveforms based on their peak and average them in groups. - Model Fitting: For each waveform (or averaged waveform): - It iterates through a range of possible surface RMS height values (`sigmah`). - It generates a theoretical nadir altimetry waveform model (`IRh1`) based on the current `sigmah` and system parameters (defined in the script, e.g., bandwidth, pulse width, wavelength, altitude). - It compares the measured waveform (`IMP`) to the theoretical model within a specified signal threshold (`cutdB`). - It calculates a metric (likely related to Maximum Likelihood Estimation, `mle1`, and Mean Integral Relative Error, `mire`) to find the best-fitting `sigmah`. - Bias Correction: Applies a hardcoded bias correction (`shift_nadir`) to the best-fit `sigmah` value. - Output: - Stores the final estimated RMS height (`SIGMAHSTIMATO`) and the fit error (`ERRORE`) for each waveform. - Generates plots comparing the measured waveform to the best-fit theoretical model for visualization. - Calculates and displays the mean and standard deviation of the estimated roughness across all processed waveforms. Input Data Files: - T23_WS.mat, T30_WS.mat, T39_WS.mat, T50_WS.mat, T56_WS.mat, T104_WS.mat, T126_WS.mat: These `.mat` files contain the pre-processed altimetry waveform data and necessary metadata required by the `NADIR_ROUGHNESS_EST_PLAINS.m` script for each respective Cassini flyby altimetry track section. Usage: - Run the `NADIR_ROUGHNESS_EST_PLAINS.m` script in MATLAB. - When prompted, select one of the `TXX_WS.mat` data files. - The script will process the waveforms in that file and display results and plots.