InSAR Product Archive

#!/usr/bin/env python3
"""
Create Multi-Product HDF5 File (InSAR Format v2.0)

This script demonstrates how to create an HDF5 file containing multiple
InSAR product types (interferograms, time series, velocity) following the
EarthScope/UNAVCO format specification version 2.0.

Replace synthetic data generation with your actual processed InSAR data.
"""

import h5py
import numpy as np
from datetime import datetime

# ============================================================================
# FILE SETUP
# ============================================================================
filename = 'insar_products.h5'
print(f"Creating HDF5 file: {filename}")

# Open HDF5 file in write mode
with h5py.File(filename, 'w') as f:
    
    # ========================================================================
    # SECTION 1: ROOT-LEVEL METADATA (applies to entire file)
    # ========================================================================
    print("\n1. Writing root-level metadata...")
    f.attrs['description'] = """InSAR results of the Los Angeles Basin covering 2020-2024. 
    It includes raw interferograms, time series and velocity products from
    ascending and descending Sentinel-1 and ALOS-2 data."""
    
    # REQUIRED: Software used for processing
    f.attrs['processing_software'] = 'ISCE2 v2.6.3 + MintPy v1.5.1'
    
    # REQUIRED: File creation timestamp (ISO 8601 format)
    f.attrs['history'] = datetime.now().isoformat()
    
    # REQUIRED: Sign convention for phase and displacement
    f.attrs['sign_convention'] = ('Negative phase change and Positive LOS '
                                   'displacement corresponds to surface motion '
                                   'toward the sensor')
    
    # RECOMMENDED: Author information as JSON string
    f.attrs['creators'] = ('[{"name":"Jane Smith","institution":"University of Example"},'
                           '{"name":"John Doe","institution":"Research Institute"}]')
    
    # RECOMMENDED: Related publication
    f.attrs['publication'] = 'Smith et al. (2024). doi:10.xxxx/yyyy'
    
    print("   ✓ Root metadata written")
    
    # ========================================================================
    # SECTION 2: CREATE FIRST TRACK (2D Grid - All Product Types)
    # ========================================================================
    print("\n2. Creating Track 1: ALOS2_073_A (2D grid with all products)...")
    
    # Create track group (naming convention: PLATFORM_ORBIT_DIRECTION)
    track1 = f.create_group('ALOS2_073_A')
    
    # ------------------------------------------------------------------------
    # 2.1: Track-Level REQUIRED Metadata
    # ------------------------------------------------------------------------
    
    # NEW in v2.0: Declare which product types this track contains
    track1.attrs['product_types'] = '["INTERFEROGRAM", "TIMESERIES", "VELOCITY"]'
    
    # NEW in v2.0: Coordinate reference system (must be EPSG:4326)
    track1.attrs['coordinate_reference_system'] = 'EPSG:4326'
    
    # Platform and acquisition parameters
    track1.attrs['platform'] = 'ALOS-2'
    track1.attrs['relative_orbit'] = 73
    track1.attrs['flight_direction'] = 'A'  # A=Ascending, D=Descending
    track1.attrs['look_direction'] = 'R'    # R=Right, L=Left
    track1.attrs['beam_mode'] = 'WD1'
    track1.attrs['wavelength'] = 0.236  # L-band in meters
    
    # Geographic extent as WKT polygon (lon, lat)
    track1.attrs['scene_footprint'] = ('POLYGON((-118.5 34.0, -118.0 34.0, '
                                        '-118.0 34.5, -118.5 34.5, -118.5 34.0))')
    
    # Temporal coverage
    track1.attrs['first_date'] = '2024-01-01'
    track1.attrs['last_date'] = '2024-04-01'
    track1.attrs['time_acquisition'] = '10:23'  # UTC time
    
    # RECOMMENDED metadata
    track1.attrs['beam_swath'] = 'W1'
    track1.attrs['polarization'] = 'VV'
    track1.attrs['atmos_correct_method'] = 'GACOS'
    track1.attrs['processing_dem'] = 'SRTM30'
    
    print("   ✓ Track metadata written")
    
    # ------------------------------------------------------------------------
    # 2.2: Geographic Coordinates (REQUIRED - stored at track level)
    # ------------------------------------------------------------------------
    print("   Creating geographic coordinates...")
    
    # Define grid dimensions
    # TODO: Replace with your actual data dimensions
    nrows, ncols = 1000, 1200
    
    # Generate coordinate arrays
    # TODO: Replace with your actual longitude/latitude arrays
    # For regular grids, create from 1D vectors:
    lon_1d = np.linspace(-118.5, -118.0, ncols, dtype=np.float32)
    lat_1d = np.linspace(34.0, 34.5, nrows, dtype=np.float32)
    lon_2d, lat_2d = np.meshgrid(lon_1d, lat_1d)
    
    # Store longitude with compression
    lon_dset = track1.create_dataset('longitude', 
                                      data=lon_2d,
                                      compression='gzip',
                                      compression_opts=6)  # 1-9, higher=smaller/slower
    lon_dset.attrs['description'] = 'Longitude coordinate'
    lon_dset.attrs['units'] = 'degrees_east'
    lon_dset.attrs['valid_range'] = np.array([-180.0, 180.0], dtype=np.float32)
    lon_dset.attrs['standard_name'] = 'longitude'
    
    # Store latitude with compression
    lat_dset = track1.create_dataset('latitude',
                                      data=lat_2d,
                                      compression='gzip',
                                      compression_opts=6)
    lat_dset.attrs['description'] = 'Latitude coordinate'
    lat_dset.attrs['units'] = 'degrees_north'
    lat_dset.attrs['valid_range'] = np.array([-90.0, 90.0], dtype=np.float32)
    lat_dset.attrs['standard_name'] = 'latitude'
    
    print(f"   ✓ Coordinates: longitude {lon_2d.shape}, latitude {lat_2d.shape}")
    
    # ------------------------------------------------------------------------
    # 2.3: Line-of-Sight (LOS) Vectors (REQUIRED - shared across products)
    # ------------------------------------------------------------------------
    print("   Creating LOS vectors...")
    
    # TODO: Replace with your actual LOS vectors from geometry calculation
    # Generate synthetic LOS vectors (East, North, Up components)
    los_e = np.random.uniform(0.35, 0.45, (nrows, ncols)).astype(np.float32)
    los_n = np.random.uniform(-0.05, 0.05, (nrows, ncols)).astype(np.float32)
    los_u = np.random.uniform(0.75, 0.85, (nrows, ncols)).astype(np.float32)
    
    # Normalize to unit vectors (important!)
    magnitude = np.sqrt(los_e**2 + los_n**2 + los_u**2)
    los_e /= magnitude
    los_n /= magnitude
    los_u /= magnitude
    
    # Store East component
    los_e_dset = track1.create_dataset('line_of_sight_e',
                                        data=los_e,
                                        compression='gzip',
                                        compression_opts=6)
    los_e_dset.attrs['description'] = 'LOS unit vector - East component'
    los_e_dset.attrs['units'] = 'dimensionless'
    
    # Store North component
    los_n_dset = track1.create_dataset('line_of_sight_n',
                                        data=los_n,
                                        compression='gzip',
                                        compression_opts=6)
    los_n_dset.attrs['description'] = 'LOS unit vector - North component'
    los_n_dset.attrs['units'] = 'dimensionless'
    
    # Store Up component
    los_u_dset = track1.create_dataset('line_of_sight_u',
                                        data=los_u,
                                        compression='gzip',
                                        compression_opts=6)
    los_u_dset.attrs['description'] = 'LOS unit vector - Up component'
    los_u_dset.attrs['units'] = 'dimensionless'
    
    print("   ✓ LOS vectors created")
    
    # ------------------------------------------------------------------------
    # 2.4: INTERFEROGRAM Product Group
    # ------------------------------------------------------------------------
    print("   Creating INTERFEROGRAM group...")
    
    # Create interferogram group
    ifg_group = track1.create_group('INTERFEROGRAM')
    
    # Create first interferogram date pair (group name: REFDATE_SECDATE)
    pair1 = ifg_group.create_group('20240101_20240113')
    
    # Date pair metadata (RECOMMENDED)
    pair1.attrs['reference_date'] = '20240101'
    pair1.attrs['secondary_date'] = '20240113'
    pair1.attrs['temporal_baseline_days'] = 12
    pair1.attrs['baseline_perp'] = 45.2  # meters
    pair1.attrs['average_coherence'] = 0.72
    pair1.attrs['percent_unwrapped'] = 94.2
    
    # Generate unwrapped phase (TODO: Replace with your actual data)
    x = np.linspace(0, 4*np.pi, ncols)
    y = np.linspace(0, 4*np.pi, nrows)
    X, Y = np.meshgrid(x, y)
    unwrapped_phase = (np.sin(X/3) * np.cos(Y/3) * 2 * np.pi).astype(np.float32)
    
    # Store unwrapped interferogram
    unwrapped = pair1.create_dataset('unwrapped_interferogram',
                                      data=unwrapped_phase,
                                      compression='gzip',
                                      compression_opts=6)
    unwrapped.attrs['description'] = 'Unwrapped interferometric phase'
    unwrapped.attrs['units'] = 'radians'
    unwrapped.attrs['unwrap_method'] = 'SNAPHU'
    
    # Generate and store wrapped phase
    wrapped_phase = np.angle(np.exp(1j * unwrapped_phase)).astype(np.float32)
    wrapped = pair1.create_dataset('wrapped_interferogram',
                                    data=wrapped_phase,
                                    compression='gzip',
                                    compression_opts=6)
    wrapped.attrs['description'] = 'Wrapped interferometric phase'
    wrapped.attrs['units'] = 'radians'
    wrapped.attrs['valid_range'] = np.array([-np.pi, np.pi], dtype=np.float32)
    
    # Generate and store coherence (TODO: Replace with your actual coherence)
    coherence = np.random.uniform(0.5, 0.9, (nrows, ncols)).astype(np.float32)
    corr = pair1.create_dataset('correlation',
                                 data=coherence,
                                 compression='gzip',
                                 compression_opts=6)
    corr.attrs['description'] = 'Interferometric coherence'
    corr.attrs['units'] = 'dimensionless'
    corr.attrs['valid_range'] = np.array([0.0, 1.0], dtype=np.float32)
    
    # Create second interferogram pair
    pair2 = ifg_group.create_group('20240113_20240125')
    pair2.attrs['reference_date'] = '20240113'
    pair2.attrs['secondary_date'] = '20240125'
    pair2.attrs['temporal_baseline_days'] = 12
    pair2.attrs['baseline_perp'] = 38.7
    
    # Store second pair data (using modified first pair as example)
    unwrapped2 = pair2.create_dataset('unwrapped_interferogram',
                                       data=(unwrapped_phase * 0.8).astype(np.float32),
                                       compression='gzip',
                                       compression_opts=6)
    unwrapped2.attrs['description'] = 'Unwrapped interferometric phase'
    unwrapped2.attrs['units'] = 'radians'
    
    print("   ✓ INTERFEROGRAM group with 2 date pairs")
    
    # ------------------------------------------------------------------------
    # 2.5: TIMESERIES Product Group
    # ------------------------------------------------------------------------
    print("   Creating TIMESERIES group...")
    
    # Create time series group
    ts_group = track1.create_group('TIMESERIES')
    
    # Group-level metadata (REQUIRED for time series)
    ts_group.attrs['reference_date'] = '20240101'
    ts_group.attrs['num_dates'] = 4
    ts_group.attrs['estimation_method'] = 'SBAS'
    ts_group.attrs['temporal_coherence_mean'] = 0.85
    
    # Create time series datasets for each date
    dates = ['20240101', '20240113', '20240125', '20240208']
    
    for idx, date_str in enumerate(dates):
        # TODO: Replace with your actual displacement data
        if idx == 0:
            # Reference date: all zeros
            displacement = np.zeros((nrows, ncols), dtype=np.float32)
        else:
            # Cumulative displacement in meters relative to reference date
            displacement = (np.random.randn(nrows, ncols) * 0.01 * idx).astype(np.float32)
        
        # Store displacement dataset (naming: dLOS_YYYYMMDD)
        dset = ts_group.create_dataset(f'dLOS_{date_str}',
                                        data=displacement,
                                        compression='gzip',
                                        compression_opts=6)
        dset.attrs['description'] = 'Cumulative LOS displacement relative to reference date'
        dset.attrs['units'] = 'meters'
        dset.attrs['acquisition_date'] = date_str
        dset.attrs['reference_date'] = '20240101'
        dset.attrs['temporal_coherence'] = 0.85 + np.random.uniform(-0.05, 0.05)
    
    print("   ✓ TIMESERIES group with 4 dates")
    
    # ------------------------------------------------------------------------
    # 2.6: VELOCITY Product Group
    # ------------------------------------------------------------------------
    print("   Creating VELOCITY group...")
    
    # Create velocity group
    vel_group = track1.create_group('VELOCITY')
    
    # Group-level metadata (RECOMMENDED)
    vel_group.attrs['time_span_start'] = '2024-01-01'
    vel_group.attrs['time_span_end'] = '2024-04-01'
    vel_group.attrs['estimation_method'] = 'linear regression'
    
    # TODO: Replace with your actual velocity data
    velocity = (np.random.randn(nrows, ncols) * 0.05).astype(np.float32)
    
    # Store velocity
    vel_dset = vel_group.create_dataset('velocity',
                                         data=velocity,
                                         compression='gzip',
                                         compression_opts=6)
    vel_dset.attrs['description'] = 'Mean LOS velocity'
    vel_dset.attrs['units'] = 'm/year'
    
    # Store velocity uncertainty (optional but recommended)
    velocity_std = (np.abs(np.random.randn(nrows, ncols)) * 0.01).astype(np.float32)
    vel_std = vel_group.create_dataset('velocity_std',
                                        data=velocity_std,
                                        compression='gzip',
                                        compression_opts=6)
    vel_std.attrs['description'] = 'Standard deviation of LOS velocity'
    vel_std.attrs['units'] = 'm/year'
    
    print("   ✓ VELOCITY group created")
    print(f"   ✓ Track 1 complete: {lon_2d.shape[0]}×{lon_2d.shape[1]} grid")
    
    # ========================================================================
    # SECTION 3: CREATE SECOND TRACK (Sparse Points - PS-InSAR Example)
    # ========================================================================
    print("\n3. Creating Track 2: S1_064_D (sparse points with velocity only)...")
    
    # Create second track
    track2 = f.create_group('S1_064_D')
    
    # This track only contains velocity (PS-InSAR results)
    track2.attrs['product_types'] = '["VELOCITY"]'
    track2.attrs['coordinate_reference_system'] = 'EPSG:4326'
    
    # Track metadata
    track2.attrs['platform'] = 'SENTINEL-1'
    track2.attrs['relative_orbit'] = 64
    track2.attrs['flight_direction'] = 'D'
    track2.attrs['look_direction'] = 'R'
    track2.attrs['beam_mode'] = 'IW'
    track2.attrs['beam_swath'] = 'IW2'
    track2.attrs['wavelength'] = 0.0555462  # C-band in meters
    track2.attrs['scene_footprint'] = ('POLYGON((-118.3 33.8, -117.8 33.8, '
                                        '-117.8 34.3, -118.3 34.3, -118.3 33.8))')
    track2.attrs['first_date'] = '2024-02-01'
    track2.attrs['last_date'] = '2024-05-15'
    track2.attrs['time_acquisition'] = '18:15'
    track2.attrs['polarization'] = 'VV'
    track2.attrs['post_processing_method'] = 'StaMPS'
    
    # ------------------------------------------------------------------------
    # 3.1: Sparse Point Coordinates (1D arrays)
    # ------------------------------------------------------------------------
    print("   Creating sparse point coordinates...")
    
    # Number of persistent scatterer points
    n_points = 15000
    
    # TODO: Replace with your actual PS point locations
    lon_points = np.random.uniform(-118.3, -117.8, n_points).astype(np.float32)
    lat_points = np.random.uniform(33.8, 34.3, n_points).astype(np.float32)
    
    # Store 1D coordinate arrays
    lon_sparse = track2.create_dataset('longitude',
                                        data=lon_points,
                                        compression='gzip',
                                        compression_opts=6)
    lon_sparse.attrs['description'] = 'Longitude coordinate'
    lon_sparse.attrs['units'] = 'degrees_east'
    lon_sparse.attrs['valid_range'] = np.array([-180.0, 180.0], dtype=np.float32)
    
    lat_sparse = track2.create_dataset('latitude',
                                        data=lat_points,
                                        compression='gzip',
                                        compression_opts=6)
    lat_sparse.attrs['description'] = 'Latitude coordinate'
    lat_sparse.attrs['units'] = 'degrees_north'
    lat_sparse.attrs['valid_range'] = np.array([-90.0, 90.0], dtype=np.float32)
    
    print(f"   ✓ Sparse coordinates: {n_points} points")
    
    # ------------------------------------------------------------------------
    # 3.2: LOS Vectors for Sparse Points (1D arrays)
    # ------------------------------------------------------------------------
    
    # TODO: Replace with your actual LOS vectors for each point
    los_e2 = np.random.uniform(-0.45, -0.35, n_points).astype(np.float32)
    los_n2 = np.random.uniform(-0.05, 0.05, n_points).astype(np.float32)
    los_u2 = np.random.uniform(0.75, 0.85, n_points).astype(np.float32)
    
    # Normalize to unit vectors
    magnitude2 = np.sqrt(los_e2**2 + los_n2**2 + los_u2**2)
    los_e2 /= magnitude2
    los_n2 /= magnitude2
    los_u2 /= magnitude2
    
    # Store LOS vectors
    los_e2_dset = track2.create_dataset('line_of_sight_e',
                                         data=los_e2,
                                         compression='gzip',
                                         compression_opts=6)
    los_e2_dset.attrs['description'] = 'LOS unit vector - East component'
    los_e2_dset.attrs['units'] = 'dimensionless'
    
    los_n2_dset = track2.create_dataset('line_of_sight_n',
                                         data=los_n2,
                                         compression='gzip',
                                         compression_opts=6)
    los_n2_dset.attrs['description'] = 'LOS unit vector - North component'
    los_n2_dset.attrs['units'] = 'dimensionless'
    
    los_u2_dset = track2.create_dataset('line_of_sight_u',
                                         data=los_u2,
                                         compression='gzip',
                                         compression_opts=6)
    los_u2_dset.attrs['description'] = 'LOS unit vector - Up component'
    los_u2_dset.attrs['units'] = 'dimensionless'
    
    # ------------------------------------------------------------------------
    # 3.3: Velocity for Sparse Points
    # ------------------------------------------------------------------------
    print("   Creating VELOCITY group for sparse points...")
    
    vel_group2 = track2.create_group('VELOCITY')
    vel_group2.attrs['time_span_start'] = '2024-02-01'
    vel_group2.attrs['time_span_end'] = '2024-05-15'
    vel_group2.attrs['estimation_method'] = 'PS-InSAR linear regression'
    
    # TODO: Replace with your actual velocity for each point
    velocity_sparse = (np.random.randn(n_points) * 0.02).astype(np.float32)
    
    vel_dset2 = vel_group2.create_dataset('velocity',
                                           data=velocity_sparse,
                                           compression='gzip',
                                           compression_opts=6)
    vel_dset2.attrs['description'] = 'Mean LOS velocity'
    vel_dset2.attrs['units'] = 'm/year'
    
    # Velocity uncertainty
    vel_std_sparse = (np.abs(np.random.randn(n_points)) * 0.005).astype(np.float32)
    vel_std2 = vel_group2.create_dataset('velocity_std',
                                          data=vel_std_sparse,
                                          compression='gzip',
                                          compression_opts=6)
    vel_std2.attrs['description'] = 'Standard deviation of LOS velocity'
    vel_std2.attrs['units'] = 'm/year'
    
    print(f"   ✓ Track 2 complete: {n_points} PS points")

# ============================================================================
# COMPLETION MESSAGE
# ============================================================================
print(f"\n{'='*70}")
print(f"✓ File created successfully: {filename}")
print(f"{'='*70}")

%% Create Multi-Product HDF5 File (InSAR Format v2.0) - MATLAB
% 
% This script demonstrates how to create an HDF5 file containing multiple
% InSAR product types (interferograms, time series, velocity) following the
% EarthScope/UNAVCO format specification version 2.0.
%
% MATLAB HDF5 API: MATLAB uses h5create, h5write, h5writeatt functions
%
% Replace synthetic data generation with your actual processed InSAR data.

clear; clc;

%% ========================================================================
%% SECTION 1: FILE SETUP
%% ========================================================================
filename = 'insar_products_matlab.h5';
fprintf('Creating HDF5 file: %s\n', filename);

% Delete if exists (to start fresh)
if exist(filename, 'file')
    delete(filename);
end

%% ========================================================================
%% SECTION 2: ROOT-LEVEL METADATA (applies to entire file)
%% ========================================================================
fprintf('\n1. Writing root-level metadata...\n');

% REQUIRED: Software used for processing
h5writeatt(filename, '/', 'processing_software', 'ISCE2 v2.6.3 + MintPy v1.5.1');

% REQUIRED: File creation timestamp (ISO 8601 format)
h5writeatt(filename, '/', 'history', datestr(now, 'yyyy-mm-ddTHH:MM:SS'));

% REQUIRED: Sign convention for phase and displacement
sign_conv = ['Negative phase change and Positive LOS displacement ' ...
             'corresponds to surface motion toward the sensor'];
h5writeatt(filename, '/', 'sign_convention', sign_conv);

% RECOMMENDED: Author information as JSON string
creators_json = '[{"name":"Jane Smith","institution":"University of Example"}]';
h5writeatt(filename, '/', 'creators', creators_json);

% RECOMMENDED: Related publication
h5writeatt(filename, '/', 'publication', 'Smith et al. (2024). doi:10.xxxx/yyyy');

% RECOMMENDED: Description of the dataset
h5writeatt(filename, '/', 'description', ...
    'InSAR results of the Los Angeles Basin covering 2020-2024. It includes raw interferograms, 
    time series and velocity products from ascending and descending Sentinel-1 and ALOS-2 data.');

fprintf('   ✓ Root metadata written\n');

%% ========================================================================
%% SECTION 3: CREATE FIRST TRACK (2D Grid - All Product Types)
%% ========================================================================
fprintf('\n2. Creating Track 1: ALOS2_073_A (2D grid with all products)...\n');

% Track name (naming convention: PLATFORM_ORBIT_DIRECTION)
track1_name = '/ALOS2_073_A';

%% ------------------------------------------------------------------------
%% 3.1: Track-Level REQUIRED Metadata
%% ------------------------------------------------------------------------

% NEW in v2.0: Declare which product types this track contains
product_types_json = '["INTERFEROGRAM", "TIMESERIES", "VELOCITY"]';
h5writeatt(filename, track1_name, 'product_types', product_types_json);

% NEW in v2.0: Coordinate reference system (must be EPSG:4326)
h5writeatt(filename, track1_name, 'coordinate_reference_system', 'EPSG:4326');

% Platform and acquisition parameters
h5writeatt(filename, track1_name, 'platform', 'ALOS-2');
h5writeatt(filename, track1_name, 'relative_orbit', int32(73));
h5writeatt(filename, track1_name, 'flight_direction', 'A'); % A=Ascending, D=Descending
h5writeatt(filename, track1_name, 'look_direction', 'R');   % R=Right, L=Left
h5writeatt(filename, track1_name, 'beam_mode', 'WD1');
h5writeatt(filename, track1_name, 'wavelength', 0.236);  % L-band in meters

% Geographic extent as WKT polygon (lon, lat)
footprint = 'POLYGON((-118.5 34.0, -118.0 34.0, -118.0 34.5, -118.5 34.5, -118.5 34.0))';
h5writeatt(filename, track1_name, 'scene_footprint', footprint);

% Temporal coverage
h5writeatt(filename, track1_name, 'first_date', '2024-01-01');
h5writeatt(filename, track1_name, 'last_date', '2024-04-01');
h5writeatt(filename, track1_name, 'time_acquisition', '10:23'); % UTC time

% RECOMMENDED metadata
h5writeatt(filename, track1_name, 'beam_swath', 'W1');
h5writeatt(filename, track1_name, 'polarization', 'VV');
h5writeatt(filename, track1_name, 'atmos_correct_method', 'GACOS');
h5writeatt(filename, track1_name, 'processing_dem', 'SRTM30');

fprintf('   ✓ Track metadata written\n');

%% ------------------------------------------------------------------------
%% 3.2: Geographic Coordinates (REQUIRED - stored at track level)
%% ------------------------------------------------------------------------
fprintf('   Creating geographic coordinates...\n');

% Define grid dimensions
% TODO: Replace with your actual data dimensions
nrows = 1000;
ncols = 1200;

% Generate coordinate arrays
% TODO: Replace with your actual longitude/latitude arrays
% For regular grids, create from 1D vectors:
lon_1d = linspace(-118.5, -118.0, ncols);
lat_1d = linspace(34.0, 34.5, nrows);
[lon_2d, lat_2d] = meshgrid(lon_1d, lat_1d);

% Convert to single precision to save space
lon_2d = single(lon_2d);
lat_2d = single(lat_2d);

% Create and write longitude dataset
% Note: MATLAB h5create requires dimensions in reverse order (MATLAB is column-major)
h5create(filename, [track1_name '/longitude'], size(lon_2d), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [track1_name '/longitude'], lon_2d);

% Write longitude attributes
h5writeatt(filename, [track1_name '/longitude'], 'description', 'Longitude coordinate');
h5writeatt(filename, [track1_name '/longitude'], 'units', 'degrees_east');
h5writeatt(filename, [track1_name '/longitude'], 'valid_range', single([-180.0, 180.0]));
h5writeatt(filename, [track1_name '/longitude'], 'standard_name', 'longitude');

% Create and write latitude dataset
h5create(filename, [track1_name '/latitude'], size(lat_2d), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [track1_name '/latitude'], lat_2d);

% Write latitude attributes
h5writeatt(filename, [track1_name '/latitude'], 'description', 'Latitude coordinate');
h5writeatt(filename, [track1_name '/latitude'], 'units', 'degrees_north');
h5writeatt(filename, [track1_name '/latitude'], 'valid_range', single([-90.0, 90.0]));
h5writeatt(filename, [track1_name '/latitude'], 'standard_name', 'latitude');

fprintf('   ✓ Coordinates: longitude [%d×%d], latitude [%d×%d]\n', ...
        size(lon_2d,1), size(lon_2d,2), size(lat_2d,1), size(lat_2d,2));

%% ------------------------------------------------------------------------
%% 3.3: Line-of-Sight (LOS) Vectors (REQUIRED - shared across products)
%% ------------------------------------------------------------------------
fprintf('   Creating LOS vectors...\n');

% TODO: Replace with your actual LOS vectors from geometry calculation
% Generate synthetic LOS vectors (East, North, Up components)
los_e = rand(nrows, ncols, 'single') * 0.1 + 0.4;  % ~0.35-0.45
los_n = rand(nrows, ncols, 'single') * 0.1 - 0.05; % ~-0.05-0.05
los_u = rand(nrows, ncols, 'single') * 0.1 + 0.8;  % ~0.75-0.85

% Normalize to unit vectors (IMPORTANT!)
magnitude = sqrt(los_e.^2 + los_n.^2 + los_u.^2);
los_e = los_e ./ magnitude;
los_n = los_n ./ magnitude;
los_u = los_u ./ magnitude;

% Create and write East component
h5create(filename, [track1_name '/line_of_sight_e'], size(los_e), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [track1_name '/line_of_sight_e'], los_e);
h5writeatt(filename, [track1_name '/line_of_sight_e'], ...
           'description', 'LOS unit vector - East component');
h5writeatt(filename, [track1_name '/line_of_sight_e'], 'units', 'dimensionless');

% Create and write North component
h5create(filename, [track1_name '/line_of_sight_n'], size(los_n), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [track1_name '/line_of_sight_n'], los_n);
h5writeatt(filename, [track1_name '/line_of_sight_n'], ...
           'description', 'LOS unit vector - North component');
h5writeatt(filename, [track1_name '/line_of_sight_n'], 'units', 'dimensionless');

% Create and write Up component
h5create(filename, [track1_name '/line_of_sight_u'], size(los_u), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [track1_name '/line_of_sight_u'], los_u);
h5writeatt(filename, [track1_name '/line_of_sight_u'], ...
           'description', 'LOS unit vector - Up component');
h5writeatt(filename, [track1_name '/line_of_sight_u'], 'units', 'dimensionless');

fprintf('   ✓ LOS vectors created\n');

%% ------------------------------------------------------------------------
%% 3.4: INTERFEROGRAM Product Group
%% ------------------------------------------------------------------------
fprintf('   Creating INTERFEROGRAM group...\n');

% Create first interferogram date pair (group name: REFDATE_SECDATE)
pair1_name = [track1_name '/INTERFEROGRAM/20240101_20240113'];

% Date pair metadata (RECOMMENDED)
h5writeatt(filename, pair1_name, 'reference_date', '20240101');
h5writeatt(filename, pair1_name, 'secondary_date', '20240113');
h5writeatt(filename, pair1_name, 'temporal_baseline_days', int32(12));
h5writeatt(filename, pair1_name, 'baseline_perp', 45.2);  % meters
h5writeatt(filename, pair1_name, 'average_coherence', 0.72);
h5writeatt(filename, pair1_name, 'percent_unwrapped', 94.2);

% Generate unwrapped phase (TODO: Replace with your actual data)
[X, Y] = meshgrid(linspace(0, 4*pi, ncols), linspace(0, 4*pi, nrows));
unwrapped_phase = single(sin(X/3) .* cos(Y/3) * 2 * pi);

% Create and write unwrapped interferogram
h5create(filename, [pair1_name '/unwrapped_interferogram'], size(unwrapped_phase), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [pair1_name '/unwrapped_interferogram'], unwrapped_phase);
h5writeatt(filename, [pair1_name '/unwrapped_interferogram'], ...
           'description', 'Unwrapped interferometric phase');
h5writeatt(filename, [pair1_name '/unwrapped_interferogram'], 'units', 'radians');
h5writeatt(filename, [pair1_name '/unwrapped_interferogram'], 'unwrap_method', 'SNAPHU');

% Generate and store wrapped phase
wrapped_phase = single(angle(exp(1i * unwrapped_phase)));

h5create(filename, [pair1_name '/wrapped_interferogram'], size(wrapped_phase), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [pair1_name '/wrapped_interferogram'], wrapped_phase);
h5writeatt(filename, [pair1_name '/wrapped_interferogram'], ...
           'description', 'Wrapped interferometric phase');
h5writeatt(filename, [pair1_name '/wrapped_interferogram'], 'units', 'radians');
h5writeatt(filename, [pair1_name '/wrapped_interferogram'], ...
           'valid_range', single([-pi, pi]));

% Generate and store coherence (TODO: Replace with your actual coherence)
coherence = single(rand(nrows, ncols) * 0.4 + 0.5);  % 0.5-0.9

h5create(filename, [pair1_name '/correlation'], size(coherence), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [pair1_name '/correlation'], coherence);
h5writeatt(filename, [pair1_name '/correlation'], ...
           'description', 'Interferometric coherence');
h5writeatt(filename, [pair1_name '/correlation'], 'units', 'dimensionless');
h5writeatt(filename, [pair1_name '/correlation'], 'valid_range', single([0.0, 1.0]));

% Create second interferogram pair
pair2_name = [track1_name '/INTERFEROGRAM/20240113_20240125'];
h5writeatt(filename, pair2_name, 'reference_date', '20240113');
h5writeatt(filename, pair2_name, 'secondary_date', '20240125');
h5writeatt(filename, pair2_name, 'temporal_baseline_days', int32(12));
h5writeatt(filename, pair2_name, 'baseline_perp', 38.7);

% Store second pair data (using modified first pair as example)
h5create(filename, [pair2_name '/unwrapped_interferogram'], size(unwrapped_phase), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, [pair2_name '/unwrapped_interferogram'], unwrapped_phase * 0.8);
h5writeatt(filename, [pair2_name '/unwrapped_interferogram'], ...
           'description', 'Unwrapped interferometric phase');
h5writeatt(filename, [pair2_name '/unwrapped_interferogram'], 'units', 'radians');

fprintf('   ✓ INTERFEROGRAM group with 2 date pairs\n');

%% ------------------------------------------------------------------------
%% 3.5: TIMESERIES Product Group
%% ------------------------------------------------------------------------
fprintf('   Creating TIMESERIES group...\n');

% Time series group path
ts_group = [track1_name '/TIMESERIES'];

% Group-level metadata (REQUIRED for time series)
h5writeatt(filename, ts_group, 'reference_date', '20240101');
h5writeatt(filename, ts_group, 'num_dates', int32(4));
h5writeatt(filename, ts_group, 'estimation_method', 'SBAS');
h5writeatt(filename, ts_group, 'temporal_coherence_mean', 0.85);

% Create time series datasets for each date
dates = {'20240101', '20240113', '20240125', '20240208'};

for idx = 1:length(dates)
    date_str = dates{idx};
    
    % TODO: Replace with your actual displacement data
    if idx == 1
        % Reference date: all zeros
        displacement = zeros(nrows, ncols, 'single');
    else
        % Cumulative displacement in meters relative to reference date
        displacement = single(randn(nrows, ncols) * 0.01 * idx);
    end
    
    % Store displacement dataset (naming: dLOS_YYYYMMDD)
    dset_path = [ts_group '/dLOS_' date_str];
    h5create(filename, dset_path, size(displacement), ...
             'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
    h5write(filename, dset_path, displacement);
    
    % Write attributes
    h5writeatt(filename, dset_path, 'description', ...
               'Cumulative LOS displacement relative to reference date');
    h5writeatt(filename, dset_path, 'units', 'meters');
    h5writeatt(filename, dset_path, 'acquisition_date', date_str);
    h5writeatt(filename, dset_path, 'reference_date', '20240101');
    h5writeatt(filename, dset_path, 'temporal_coherence', ...
               0.85 + (rand - 0.5) * 0.1);
end

fprintf('   ✓ TIMESERIES group with %d dates\n', length(dates));

%% ------------------------------------------------------------------------
%% 3.6: VELOCITY Product Group
%% ------------------------------------------------------------------------
fprintf('   Creating VELOCITY group...\n');

% Velocity group path
vel_group = [track1_name '/VELOCITY'];

% Group-level metadata (RECOMMENDED)
h5writeatt(filename, vel_group, 'time_span_start', '2024-01-01');
h5writeatt(filename, vel_group, 'time_span_end', '2024-04-01');
h5writeatt(filename, vel_group, 'estimation_method', 'linear regression');

% TODO: Replace with your actual velocity data
velocity = single(randn(nrows, ncols) * 0.05);

% Store velocity
vel_path = [vel_group '/velocity'];
h5create(filename, vel_path, size(velocity), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, vel_path, velocity);
h5writeatt(filename, vel_path, 'description', 'Mean LOS velocity');
h5writeatt(filename, vel_path, 'units', 'm/year');

% Store velocity uncertainty (optional but recommended)
velocity_std = single(abs(randn(nrows, ncols)) * 0.01);
vel_std_path = [vel_group '/velocity_std'];
h5create(filename, vel_std_path, size(velocity_std), ...
         'Datatype', 'single', 'ChunkSize', [100 100], 'Deflate', 6);
h5write(filename, vel_std_path, velocity_std);
h5writeatt(filename, vel_std_path, 'description', ...
           'Standard deviation of LOS velocity');
h5writeatt(filename, vel_std_path, 'units', 'm/year');

fprintf('   ✓ VELOCITY group created\n');
fprintf('   ✓ Track 1 complete: %d×%d grid\n', nrows, ncols);

%% ========================================================================
%% SECTION 4: COMPLETION MESSAGE
%% ========================================================================
fprintf('\n%s\n', repmat('=', 1, 70));
fprintf('✓ File created successfully: %s\n', filename);
fprintf('%s\n\n', repmat('=', 1, 70));
fprintf('File contents:\n');
fprintf('  Track 1 (ALOS2_073_A): 2D grid (%d×%d)\n', nrows, ncols);
fprintf('    - INTERFEROGRAM: 2 date pairs\n');
fprintf('    - TIMESERIES: %d dates\n', length(dates));
fprintf('    - VELOCITY: mean velocity + uncertainty\n');
fprintf('\nExplore with: h5disp(''%s'')\n', filename);
fprintf('\nNOTE: This file contains SYNTHETIC DEMO DATA.\n');
fprintf('Replace data generation sections with your actual processed InSAR data.\n');

%% Read Multi-Product HDF5 File (InSAR Format v2.0) - MATLAB
%% Read Multi-Product HDF5 File (InSAR Format v2.0) - MATLAB
%
% This script demonstrates how to read InSAR products from HDF5 files
% following the EarthScope/UNAVCO format specification version 2.0.
%
% It shows how to:
% - Read file and track metadata
% - Extract geographic coordinates
% - Access interferograms, time series, and velocity products
% - Handle both 2D grids and sparse point data
%
% MATLAB HDF5 API: Uses h5read, h5readatt, h5info functions

clear; clc;

%% ========================================================================
%% SECTION 1: FILE SETUP
%% ========================================================================
filename = 'insar_products_matlab.h5';
fprintf('Reading HDF5 file: %s\n\n', filename);

% Check if file exists
if ~exist(filename, 'file')
    error('File not found: %s', filename);
end

%% ========================================================================
%% SECTION 2: READ ROOT-LEVEL METADATA
%% ========================================================================
fprintf('%s\n', repmat('=', 1, 70));
fprintf('FILE-WIDE METADATA\n');
fprintf('%s\n', repmat('=', 1, 70));

% Read required metadata
software = h5readatt(filename, '/', 'processing_software');
history = h5readatt(filename, '/', 'history');
sign_convention = h5readatt(filename, '/', 'sign_convention');

fprintf('Processing Software: %s\n', software);
fprintf('File Created: %s\n', history);
fprintf('Sign Convention: %s\n', sign_convention);

% Read optional metadata if present
try
    creators = h5readatt(filename, '/', 'creators');
    fprintf('\nCreators: %s\n', creators);
catch
    % Attribute doesn't exist
end

try
    publication = h5readatt(filename, '/', 'publication');
    fprintf('Publication: %s\n', publication);
catch
    % Attribute doesn't exist
end

%% ========================================================================
%% SECTION 3: GET TRACK INFORMATION
%% ========================================================================

% Get file structure
fileinfo = h5info(filename);

% Extract track names (top-level groups)
track_names = {fileinfo.Groups.Name};
n_tracks = length(track_names);

fprintf('\nFound %d track(s):\n', n_tracks);
for i = 1:n_tracks
    % Remove leading '/' from name
    track_display = track_names{i}(2:end);
    fprintf('  - %s\n', track_display);
end

%% ========================================================================
%% SECTION 4: PROCESS EACH TRACK
%% ========================================================================

for t = 1:n_tracks
    track_path = track_names{t};
    track_display = track_path(2:end);  % Remove leading '/'
    
    fprintf('\n%s\n', repmat('=', 1, 70));
    fprintf('TRACK: %s\n', track_display);
    fprintf('%s\n', repmat('=', 1, 70));
    
    %% --------------------------------------------------------------------
    %% 4.1: Read Track-Level Metadata
    %% --------------------------------------------------------------------
    fprintf('\nTrack Metadata:\n');
    
    % Parse product types (JSON array - need to parse manually in MATLAB)
    product_types_json = h5readatt(filename, track_path, 'product_types');
    fprintf('  Product Types (JSON): %s\n', product_types_json);
    
    % Simple JSON parsing for product types (extract values between quotes)
    % For more complex JSON, use jsondecode() in MATLAB R2016b+
    has_ifg = contains(product_types_json, 'INTERFEROGRAM');
    has_ts = contains(product_types_json, 'TIMESERIES');
    has_vel = contains(product_types_json, 'VELOCITY');
    
    fprintf('  Product Types: ');
    if has_ifg, fprintf('INTERFEROGRAM '); end
    if has_ts, fprintf('TIMESERIES '); end
    if has_vel, fprintf('VELOCITY'); end
    fprintf('\n');
    
    % Coordinate reference system (NEW in v2.0)
    crs = h5readatt(filename, track_path, 'coordinate_reference_system');
    fprintf('  CRS: %s\n', crs);
    
    % Platform and acquisition information
    platform = h5readatt(filename, track_path, 'platform');
    orbit = h5readatt(filename, track_path, 'relative_orbit');
    direction = h5readatt(filename, track_path, 'flight_direction');
    beam_mode = h5readatt(filename, track_path, 'beam_mode');
    wavelength = h5readatt(filename, track_path, 'wavelength');
    first_date = h5readatt(filename, track_path, 'first_date');
    last_date = h5readatt(filename, track_path, 'last_date');
    acq_time = h5readatt(filename, track_path, 'time_acquisition');
    
    fprintf('  Platform: %s\n', platform);
    fprintf('  Orbit: %d\n', orbit);
    fprintf('  Direction: %s\n', direction);
    fprintf('  Beam Mode: %s\n', beam_mode);
    fprintf('  Wavelength: %.6f m\n', wavelength);
    fprintf('  Time Span: %s to %s\n', first_date, last_date);
    fprintf('  Acquisition Time: %s UTC\n', acq_time);
    
    %% --------------------------------------------------------------------
    %% 4.2: Read Geographic Coordinates
    %% --------------------------------------------------------------------
    fprintf('\nGeographic Coordinates:\n');
    
    % Read longitude and latitude arrays
    longitude = h5read(filename, [track_path '/longitude']);
    latitude = h5read(filename, [track_path '/latitude']);
    
    fprintf('  Longitude: size=[%s], range=[%.4f, %.4f]°\n', ...
            mat2str(size(longitude)), min(longitude(:)), max(longitude(:)));
    fprintf('  Latitude: size=[%s], range=[%.4f, %.4f]°\n', ...
            mat2str(size(latitude)), min(latitude(:)), max(latitude(:)));
    
    % Determine coordinate type
    if ndims(longitude) == 2 && all(size(longitude) > 1)
        coord_type = '2D regular grid';
        n_pixels = numel(longitude);
    elseif isvector(longitude)
        coord_type = '1D sparse points';
        n_pixels = length(longitude);
    else
        coord_type = 'Unknown dimension';
        n_pixels = numel(longitude);
    end
    
    fprintf('  Type: %s (%s points/pixels)\n', coord_type, ...
            num2str(n_pixels, '%d'));
    
    % Read LOS vectors (shared across all products)
    los_e = h5read(filename, [track_path '/line_of_sight_e']);
    los_n = h5read(filename, [track_path '/line_of_sight_n']);
    los_u = h5read(filename, [track_path '/line_of_sight_u']);
    
    fprintf('  LOS Vectors: size=[%s]\n', mat2str(size(los_e)));
    
    % Verify LOS vectors are unit vectors
    los_magnitude = sqrt(los_e.^2 + los_n.^2 + los_u.^2);
    fprintf('  LOS Magnitude Check: min=%.6f, max=%.6f (should be ~1.0)\n', ...
            min(los_magnitude(:)), max(los_magnitude(:)));
    
    %% --------------------------------------------------------------------
    %% 4.3: Read INTERFEROGRAM Products (if present)
    %% --------------------------------------------------------------------
    if has_ifg
        fprintf('\n--- INTERFEROGRAM PRODUCTS ---\n');
        
        % Get interferogram group info
        try
            ifg_path = [track_path '/INTERFEROGRAM'];
            ifg_info = h5info(filename, ifg_path);
            date_pairs = {ifg_info.Groups.Name};
            
            fprintf('Found %d interferogram(s)\n', length(date_pairs));
            
            % Read each date pair
            for p = 1:length(date_pairs)
                pair_path = date_pairs{p};
                % Extract date pair name (last part of path)
                parts = strsplit(pair_path, '/');
                pair_name = parts{end};
                
                fprintf('\n  Date Pair: %s\n', pair_name);
                
                % Read metadata if available
                try
                    ref_date = h5readatt(filename, pair_path, 'reference_date');
                    sec_date = h5readatt(filename, pair_path, 'secondary_date');
                    fprintf('    Reference: %s, Secondary: %s\n', ref_date, sec_date);
                catch
                    % Metadata not available
                end
                
                try
                    temp_baseline = h5readatt(filename, pair_path, 'temporal_baseline_days');
                    fprintf('    Temporal Baseline: %d days\n', temp_baseline);
                catch
                end
                
                try
                    perp_baseline = h5readatt(filename, pair_path, 'baseline_perp');
                    fprintf('    Perpendicular Baseline: %.2f m\n', perp_baseline);
                catch
                end
                
                % Read unwrapped phase
                try
                    unwrapped = h5read(filename, [pair_path '/unwrapped_interferogram']);
                    fprintf('    Unwrapped Phase: size=[%s], range=[%.2f, %.2f] rad\n', ...
                            mat2str(size(unwrapped)), min(unwrapped(:)), max(unwrapped(:)));
                    
                    % Verify dimensions match coordinates
                    if isequal(size(unwrapped), size(longitude))
                        fprintf('    ✓ Dimensions match coordinates\n');
                    else
                        fprintf('    ⚠ WARNING: Dimensions mismatch!\n');
                    end
                catch ME
                    fprintf('    Could not read unwrapped phase: %s\n', ME.message);
                end
                
                % Read wrapped phase
                try
                    wrapped = h5read(filename, [pair_path '/wrapped_interferogram']);
                    fprintf('    Wrapped Phase: size=[%s], range=[%.2f, %.2f] rad\n', ...
                            mat2str(size(wrapped)), min(wrapped(:)), max(wrapped(:)));
                catch
                    % Dataset not available
                end
                
                % Read coherence
                try
                    correlation = h5read(filename, [pair_path '/correlation']);
                    fprintf('    Coherence: mean=%.3f, std=%.3f\n', ...
                            mean(correlation(:)), std(correlation(:)));
                    
                    % Count high-coherence pixels
                    high_coh = sum(correlation(:) > 0.7);
                    pct_high_coh = 100.0 * high_coh / numel(correlation);
                    fprintf('    High Coherence (>0.7): %.1f%%\n', pct_high_coh);
                catch
                    % Dataset not available
                end
            end
        catch ME
            fprintf('  ⚠ WARNING: INTERFEROGRAM declared but could not read: %s\n', ...
                    ME.message);
        end
    end
    
    %% --------------------------------------------------------------------
    %% 4.4: Read TIMESERIES Products (if present)
    %% --------------------------------------------------------------------
    if has_ts
        fprintf('\n--- TIME SERIES PRODUCTS ---\n');
        
        try
            ts_path = [track_path '/TIMESERIES'];
            
            % Read group-level metadata
            ref_date = h5readatt(filename, ts_path, 'reference_date');
            fprintf('Reference Date: %s\n', ref_date);
            
            try
                method = h5readatt(filename, ts_path, 'estimation_method');
                fprintf('Estimation Method: %s\n', method);
            catch
            end
            
            try
                temp_coh = h5readatt(filename, ts_path, 'temporal_coherence_mean');
                fprintf('Mean Temporal Coherence: %.3f\n', temp_coh);
            catch
            end
            
            % Find all time series datasets
            ts_info = h5info(filename, ts_path);
            date_datasets = {ts_info.Datasets.Name};
            
            % Extract dates from dataset names
            dates = cell(length(date_datasets), 1);
            for d = 1:length(date_datasets)
                % Remove 'dLOS_' prefix
                dates{d} = strrep(date_datasets{d}, 'dLOS_', '');
            end
            dates = sort(dates);
            
            fprintf('\nFound %d date(s)\n', length(dates));
            fprintf('Date Range: %s to %s\n', dates{1}, dates{end});
            
            fprintf('\nSample Dates:\n');
            
            % Read reference date (should be zeros)
            ref_dataset = ['dLOS_' ref_date];
            try
                ref_disp = h5read(filename, [ts_path '/' ref_dataset]);
                fprintf('\n  %s (reference):\n', ref_date);
                fprintf('    Size: [%s]\n', mat2str(size(ref_disp)));
                fprintf('    Range: [%.6f, %.6f] m\n', ...
                        min(ref_disp(:)), max(ref_disp(:)));
                
                if all(abs(ref_disp(:)) < 1e-6)
                    fprintf('    ✓ Reference date is all zeros (correct)\n');
                else
                    fprintf('    ⚠ WARNING: Reference date not zeros!\n');
                end
            catch
                fprintf('  Could not read reference date\n');
            end
            
            % Read a few other dates
            n_samples = min(3, length(dates));
            for d = 2:n_samples
                date_str = dates{d};
                dataset_name = ['dLOS_' date_str];
                
                try
                    displacement = h5read(filename, [ts_path '/' dataset_name]);
                    
                    fprintf('\n  %s:\n', date_str);
                    fprintf('    Size: [%s]\n', mat2str(size(displacement)));
                    fprintf('    Displacement: mean=%.4f m, std=%.4f m\n', ...
                            mean(displacement(:)), std(displacement(:)));
                    fprintf('    Range: [%.4f, %.4f] m\n', ...
                            min(displacement(:)), max(displacement(:)));
                    
                    % Check for valid data
                    valid_pixels = sum(isfinite(displacement(:)));
                    pct_valid = 100.0 * valid_pixels / numel(displacement);
                    fprintf('    Valid Pixels: %.1f%%\n', pct_valid);
                catch
                    fprintf('  Could not read %s\n', date_str);
                end
            end
            
            % Example: Extract time series at a specific location (for 2D grids)
            if ndims(longitude) == 2
                fprintf('\nExample: Time Series at Specific Location\n');
                
                % Target location
                target_lon = -118.25;
                target_lat = 34.25;
                
                % Find nearest pixel
                lon_diff = abs(longitude - target_lon);
                lat_diff = abs(latitude - target_lat);
                distance = sqrt(lon_diff.^2 + lat_diff.^2);
                [~, idx] = min(distance(:));
                [row, col] = ind2sub(size(distance), idx);
                
                fprintf('  Target: (%.4f°, %.4f°)\n', target_lon, target_lat);
                fprintf('  Nearest Pixel: (%.4f°, %.4f°)\n', ...
                        longitude(row, col), latitude(row, col));
                fprintf('  Array Index: [%d, %d]\n', row, col);
                
                % Extract time series at this point
                fprintf('\n  Displacement Time Series:\n');
                n_show = min(5, length(dates));
                for d = 1:n_show
                    date_str = dates{d};
                    dataset_name = ['dLOS_' date_str];
                    
                    try
                        disp_data = h5read(filename, [ts_path '/' dataset_name]);
                        disp_value = disp_data(row, col);
                        fprintf('    %s: %+.4f m\n', date_str, disp_value);
                    catch
                        fprintf('    %s: (could not read)\n', date_str);
                    end
                end
                
                if length(dates) > n_show
                    fprintf('    ... (%d more dates)\n', length(dates) - n_show);
                end
            end
        catch ME
            fprintf('  ⚠ WARNING: TIMESERIES declared but could not read: %s\n', ...
                    ME.message);
        end
    end
    
    %% --------------------------------------------------------------------
    %% 4.5: Read VELOCITY Products (if present)
    %% --------------------------------------------------------------------
    if has_vel
        fprintf('\n--- VELOCITY PRODUCTS ---\n');
        
        try
            vel_path = [track_path '/VELOCITY'];
            
            % Read group-level metadata
            try
                time_start = h5readatt(filename, vel_path, 'time_span_start');
                time_end = h5readatt(filename, vel_path, 'time_span_end');
                fprintf('Time Span: %s to %s\n', time_start, time_end);
            catch
            end
            
            try
                method = h5readatt(filename, vel_path, 'estimation_method');
                fprintf('Estimation Method: %s\n', method);
            catch
            end
            
            % Read velocity
            velocity = h5read(filename, [vel_path '/velocity']);
            
            fprintf('\nVelocity Statistics:\n');
            fprintf('  Size: [%s]\n', mat2str(size(velocity)));
            fprintf('  Dimensions match coordinates: %d\n', ...
                    isequal(size(velocity), size(longitude)));
            fprintf('  Mean: %.4f m/year\n', mean(velocity(:)));
            fprintf('  Std Dev: %.4f m/year\n', std(velocity(:)));
            fprintf('  Range: [%.4f, %.4f] m/year\n', ...
                    min(velocity(:)), max(velocity(:)));
            
            % Count deforming pixels (|velocity| > threshold)
            threshold = 0.01;  % m/year
            deforming = abs(velocity) > threshold;
            n_deforming = sum(deforming(:));
            pct_deforming = 100.0 * n_deforming / numel(velocity);
            fprintf('  Pixels with |v| > %.2f m/year: %d (%.1f%%)\n', ...
                    threshold, n_deforming, pct_deforming);
            
            % Read velocity uncertainty if available
            try
                vel_std = h5read(filename, [vel_path '/velocity_std']);
                fprintf('\nVelocity Uncertainty:\n');
                fprintf('  Mean: %.4f m/year\n', mean(vel_std(:)));
                fprintf('  Range: [%.4f, %.4f] m/year\n', ...
                        min(vel_std(:)), max(vel_std(:)));
                
                % Calculate signal-to-noise ratio
                snr = abs(velocity) ./ vel_std;
                fprintf('  Mean SNR: %.2f\n', mean(snr(:), 'omitnan'));
            catch
                fprintf('\nVelocity uncertainty not available\n');
            end
            
            % Example: Find pixels with significant deformation
            if ndims(longitude) == 2
                fprintf('\nExample: Significant Deformation Pixels\n');
                
                % Define significance threshold
                try
                    vel_std = h5read(filename, [vel_path '/velocity_std']);
                    significant = abs(velocity) > (2 * vel_std);
                catch
                    significant = abs(velocity) > 0.02;  % 2 cm/year
                end
                
                n_sig = sum(significant(:));
                fprintf('  Significant pixels: %d (%.1f%%)\n', ...
                        n_sig, 100.0*n_sig/numel(velocity));
                
                if n_sig > 0
                    % Find location of max subsidence and uplift
                    [min_vel, min_idx] = min(velocity(:));
                    [max_vel, max_idx] = max(velocity(:));
                    
                    [min_row, min_col] = ind2sub(size(velocity), min_idx);
                    [max_row, max_col] = ind2sub(size(velocity), max_idx);
                    
                    fprintf('\n  Maximum Subsidence:\n');
                    fprintf('    Location: (%.4f°, %.4f°)\n', ...
                            longitude(min_row, min_col), latitude(min_row, min_col));
                    fprintf('    Velocity: %.4f m/year\n', min_vel);
                    
                    fprintf('\n  Maximum Uplift:\n');
                    fprintf('    Location: (%.4f°, %.4f°)\n', ...
                            longitude(max_row, max_col), latitude(max_row, max_col));
                    fprintf('    Velocity: %.4f m/year\n', max_vel);
                end
            end
            
            % For sparse points, show example points
            if isvector(longitude)
                fprintf('\nExample Points (showing first 5):\n');
                n_show = min(5, length(longitude));
                for i = 1:n_show
                    fprintf('  Point %d: lon=%.4f°, lat=%.4f°, velocity=%+.4f m/year\n', ...
                            i, longitude(i), latitude(i), velocity(i));
                end
            end
        catch ME
            fprintf('  ⚠ WARNING: VELOCITY declared but could not read: %s\n', ...
                    ME.message);
        end
    end
end

%% ========================================================================
%% COMPLETION MESSAGE
%% ========================================================================
fprintf('\n%s\n', repmat('=', 1, 70));
fprintf('✓ File read successfully\n');
fprintf('%s\n', repmat('=', 1, 70));
fprintf('\nTo explore the file structure, use:\n');
fprintf('  h5disp(''%s'')\n', filename);
fprintf('\nTo browse interactively, use:\n');
fprintf('  h5tool(''%s'')\n', filename);