Skip to content

I/O

Export and import functions for GeoTIFF and NetCDF formats.

GeoTIFF Export

Note

GeoTIFF functions require rasterio. Install with pip install eq-insar[geotiff].

save_geotiff 

save_geotiff(data: ndarray, filepath: Union[str, Path], bounds_km: Optional[Tuple[float, float, float, float]] = None, crs: str = 'EPSG:32633', pixel_size_km: float = 0.5, nodata: float = np.nan, description: str = '', metadata: Optional[Dict] = None) -> None

Save a 2D array to GeoTIFF format.

PARAMETER DESCRIPTION
data

2D array to save (ny, nx)

TYPE: ndarray

filepath

Output file path

TYPE: str or Path

bounds_km

(xmin, ymin, xmax, ymax) in km. If None, centers data at origin.

TYPE: tuple DEFAULT: None

crs

Coordinate reference system (default: UTM 33N)

TYPE: str DEFAULT: 'EPSG:32633'

pixel_size_km

Pixel size in km

TYPE: float DEFAULT: 0.5

nodata

NoData value

TYPE: float DEFAULT: nan

description

Band description

TYPE: str DEFAULT: ''

metadata

Additional metadata to include

TYPE: dict DEFAULT: None

RAISES DESCRIPTION
ImportError

If rasterio is not installed
Source code in src/eq_insar/io/geotiff.py 
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
def save_geotiff(
    data: np.ndarray,
    filepath: Union[str, Path],
    bounds_km: Optional[Tuple[float, float, float, float]] = None,
    crs: str = "EPSG:32633",  # UTM zone 33N as default
    pixel_size_km: float = 0.5,
    nodata: float = np.nan,
    description: str = "",
    metadata: Optional[Dict] = None,
) -> None:
    """
    Save a 2D array to GeoTIFF format.

    Parameters
    ----------
    data : np.ndarray
        2D array to save (ny, nx)
    filepath : str or Path
        Output file path
    bounds_km : tuple, optional
        (xmin, ymin, xmax, ymax) in km. If None, centers data at origin.
    crs : str
        Coordinate reference system (default: UTM 33N)
    pixel_size_km : float
        Pixel size in km
    nodata : float
        NoData value
    description : str
        Band description
    metadata : dict, optional
        Additional metadata to include

    Raises
    ------
    ImportError
        If rasterio is not installed
    """
    if not _check_rasterio():
        raise ImportError(
            "rasterio is required for GeoTIFF export. "
            "Install with: pip install rasterio"
        )

    import rasterio
    from rasterio.transform import from_bounds

    filepath = Path(filepath)
    filepath.parent.mkdir(parents=True, exist_ok=True)

    ny, nx = data.shape

    # Calculate bounds if not provided
    if bounds_km is None:
        half_x = nx * pixel_size_km / 2
        half_y = ny * pixel_size_km / 2
        bounds_km = (-half_x, -half_y, half_x, half_y)

    # Convert km to meters for geotransform
    xmin, ymin, xmax, ymax = [b * 1000 for b in bounds_km]

    # Create transform
    transform = from_bounds(xmin, ymin, xmax, ymax, nx, ny)

    # Prepare metadata
    profile = {
        "driver": "GTiff",
        "dtype": data.dtype,
        "width": nx,
        "height": ny,
        "count": 1,
        "crs": crs,
        "transform": transform,
        "nodata": nodata,
        "compress": "lzw",
    }

    # Write file
    with rasterio.open(filepath, "w", **profile) as dst:
        dst.write(data, 1)
        if description:
            dst.set_band_description(1, description)

        # Add custom metadata
        if metadata:
            dst.update_tags(**{str(k): str(v) for k, v in metadata.items()})

save_displacement_geotiff 

save_displacement_geotiff(result: Dict, output_dir: Union[str, Path], prefix: str = 'eq_insar', include_components: bool = True) -> Dict[str, Path]

Save displacement fields from EQ-INSAR result to GeoTIFF files.

PARAMETER DESCRIPTION
result

Output from generate_synthetic_insar or generate_timeseries

TYPE: dict

output_dir

Output directory

TYPE: str or Path

prefix

Filename prefix

TYPE: str DEFAULT: 'eq_insar'

include_components

If True, save Ue, Un, Uz components in addition to LOS

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
files

Dictionary mapping field names to output file paths

TYPE: dict

Source code in src/eq_insar/io/geotiff.py 
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
def save_displacement_geotiff(
    result: Dict,
    output_dir: Union[str, Path],
    prefix: str = "eq_insar",
    include_components: bool = True,
) -> Dict[str, Path]:
    """
    Save displacement fields from EQ-INSAR result to GeoTIFF files.

    Parameters
    ----------
    result : dict
        Output from generate_synthetic_insar or generate_timeseries
    output_dir : str or Path
        Output directory
    prefix : str
        Filename prefix
    include_components : bool
        If True, save Ue, Un, Uz components in addition to LOS

    Returns
    -------
    files : dict
        Dictionary mapping field names to output file paths
    """
    if not _check_rasterio():
        raise ImportError(
            "rasterio is required for GeoTIFF export. "
            "Install with: pip install rasterio"
        )

    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    files = {}
    meta = result.get("metadata", {})

    # Get bounds from coordinate grids
    X_km = result["X_km"]
    Y_km = result["Y_km"]
    bounds_km = (X_km.min(), Y_km.min(), X_km.max(), Y_km.max())
    pixel_size = meta.get("grid_spacing_km", 0.5)

    # Base metadata
    base_meta = {
        "source": "EQ-INSAR",
        "Mw": meta.get("Mw", "N/A"),
        "depth_km": meta.get("depth_km", "N/A"),
        "strike_deg": meta.get("strike_deg", "N/A"),
        "dip_deg": meta.get("dip_deg", "N/A"),
        "rake_deg": meta.get("rake_deg", "N/A"),
    }

    # Save LOS displacement
    filepath = output_dir / f"{prefix}_los_displacement_m.tif"
    save_geotiff(
        result["los_displacement"].astype(np.float32),
        filepath,
        bounds_km=bounds_km,
        pixel_size_km=pixel_size,
        description="LOS displacement (m, positive toward satellite)",
        metadata=base_meta,
    )
    files["los_displacement"] = filepath

    # Save displacement components
    if include_components:
        for comp, desc in [
            ("Ue", "East displacement (m)"),
            ("Un", "North displacement (m)"),
            ("Uz", "Vertical displacement (m, positive up)"),
        ]:
            if comp in result:
                filepath = output_dir / f"{prefix}_{comp}_m.tif"
                save_geotiff(
                    result[comp].astype(np.float32),
                    filepath,
                    bounds_km=bounds_km,
                    pixel_size_km=pixel_size,
                    description=desc,
                    metadata=base_meta,
                )
                files[comp] = filepath

    return files

save_phase_geotiff 

save_phase_geotiff(result: Dict, output_dir: Union[str, Path], prefix: str = 'eq_insar', include_all: bool = True) -> Dict[str, Path]

Save phase fields from EQ-INSAR result to GeoTIFF files.

PARAMETER DESCRIPTION
result

Output from generate_synthetic_insar

TYPE: dict

output_dir

Output directory

TYPE: str or Path

prefix

Filename prefix

TYPE: str DEFAULT: 'eq_insar'

include_all

If True, save unwrapped, wrapped, and noisy phases

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
files

Dictionary mapping field names to output file paths

TYPE: dict

Source code in src/eq_insar/io/geotiff.py 
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
def save_phase_geotiff(
    result: Dict,
    output_dir: Union[str, Path],
    prefix: str = "eq_insar",
    include_all: bool = True,
) -> Dict[str, Path]:
    """
    Save phase fields from EQ-INSAR result to GeoTIFF files.

    Parameters
    ----------
    result : dict
        Output from generate_synthetic_insar
    output_dir : str or Path
        Output directory
    prefix : str
        Filename prefix
    include_all : bool
        If True, save unwrapped, wrapped, and noisy phases

    Returns
    -------
    files : dict
        Dictionary mapping field names to output file paths
    """
    if not _check_rasterio():
        raise ImportError(
            "rasterio is required for GeoTIFF export. "
            "Install with: pip install rasterio"
        )

    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    files = {}
    meta = result.get("metadata", {})

    X_km = result["X_km"]
    Y_km = result["Y_km"]
    bounds_km = (X_km.min(), Y_km.min(), X_km.max(), Y_km.max())
    pixel_size = meta.get("grid_spacing_km", 0.5)

    base_meta = {
        "source": "EQ-INSAR",
        "wavelength_m": meta.get("wavelength_m", "N/A"),
    }

    phase_fields = [
        ("phase_wrapped", "Wrapped phase (rad)"),
    ]

    if include_all:
        phase_fields.extend([
            ("phase_unwrapped", "Unwrapped phase (rad)"),
            ("phase_noisy", "Wrapped phase with noise (rad)"),
        ])

    for field, desc in phase_fields:
        if field in result:
            filepath = output_dir / f"{prefix}_{field}.tif"
            save_geotiff(
                result[field].astype(np.float32),
                filepath,
                bounds_km=bounds_km,
                pixel_size_km=pixel_size,
                description=desc,
                metadata=base_meta,
            )
            files[field] = filepath

    return files

NetCDF Export/Load

Note

NetCDF functions require netCDF4. Install with pip install eq-insar[netcdf].

save_netcdf 

save_netcdf(result: Dict, filepath: Union[str, Path], include_metadata: bool = True, compression_level: int = 4) -> None

Save EQ-INSAR result to NetCDF format.

Creates a CF-compliant NetCDF file with all displacement and phase fields, plus metadata.

PARAMETER DESCRIPTION
result

Output from generate_synthetic_insar

TYPE: dict

filepath

Output file path

TYPE: str or Path

include_metadata

Include source parameters as global attributes

TYPE: bool DEFAULT: True

compression_level

Compression level (0-9, higher = smaller files)

TYPE: int DEFAULT: 4

RAISES DESCRIPTION
ImportError

If netCDF4 is not installed
Source code in src/eq_insar/io/netcdf.py 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def save_netcdf(
    result: Dict,
    filepath: Union[str, Path],
    include_metadata: bool = True,
    compression_level: int = 4,
) -> None:
    """
    Save EQ-INSAR result to NetCDF format.

    Creates a CF-compliant NetCDF file with all displacement and
    phase fields, plus metadata.

    Parameters
    ----------
    result : dict
        Output from generate_synthetic_insar
    filepath : str or Path
        Output file path
    include_metadata : bool
        Include source parameters as global attributes
    compression_level : int
        Compression level (0-9, higher = smaller files)

    Raises
    ------
    ImportError
        If netCDF4 is not installed
    """
    if not _check_netcdf4():
        raise ImportError(
            "netCDF4 is required for NetCDF export. "
            "Install with: pip install netCDF4"
        )

    import netCDF4

    filepath = Path(filepath)
    filepath.parent.mkdir(parents=True, exist_ok=True)

    X_km = result["X_km"]
    Y_km = result["Y_km"]
    meta = result.get("metadata", {})

    ny, nx = X_km.shape

    # Create file
    with netCDF4.Dataset(filepath, "w", format="NETCDF4") as nc:
        # Global attributes
        nc.title = "EQ-INSAR Synthetic InSAR Data"
        nc.institution = "EQ-INSAR Package"
        nc.source = "Synthetic earthquake deformation model"
        nc.history = f"Created {datetime.now().isoformat()}"
        nc.Conventions = "CF-1.8"

        if include_metadata:
            nc.earthquake_Mw = meta.get("Mw", np.nan)
            nc.earthquake_M0_Nm = meta.get("M0_Nm", np.nan)
            nc.earthquake_depth_km = meta.get("depth_km", np.nan)
            nc.earthquake_strike_deg = meta.get("strike_deg", np.nan)
            nc.earthquake_dip_deg = meta.get("dip_deg", np.nan)
            nc.earthquake_rake_deg = meta.get("rake_deg", np.nan)
            nc.earthquake_xcen_km = meta.get("xcen_km", np.nan)
            nc.earthquake_ycen_km = meta.get("ycen_km", np.nan)
            nc.source_type = meta.get("source_type", "unknown")
            nc.satellite = meta.get("satellite", "unknown")
            nc.wavelength_m = meta.get("wavelength_m", np.nan)
            nc.incidence_deg = meta.get("incidence_deg", np.nan)
            nc.heading_deg = meta.get("heading_deg", np.nan)

        # Create dimensions
        nc.createDimension("y", ny)
        nc.createDimension("x", nx)

        # Create coordinate variables
        x_var = nc.createVariable("x", "f4", ("x",))
        x_var.units = "km"
        x_var.long_name = "Easting"
        x_var.standard_name = "projection_x_coordinate"
        x_var[:] = X_km[0, :]

        y_var = nc.createVariable("y", "f4", ("y",))
        y_var.units = "km"
        y_var.long_name = "Northing"
        y_var.standard_name = "projection_y_coordinate"
        y_var[:] = Y_km[:, 0]

        # Compression settings
        comp = {"zlib": True, "complevel": compression_level}

        # Data variables
        variables = [
            ("los_displacement", "m", "LOS displacement", "positive toward satellite"),
            ("Ue", "m", "East displacement", "positive eastward"),
            ("Un", "m", "North displacement", "positive northward"),
            ("Uz", "m", "Vertical displacement", "positive upward"),
            ("phase_unwrapped", "rad", "Unwrapped interferometric phase", ""),
            ("phase_wrapped", "rad", "Wrapped interferometric phase", ""),
            ("phase_noisy", "rad", "Wrapped phase with noise", ""),
        ]

        for name, units, long_name, comment in variables:
            if name in result:
                var = nc.createVariable(name, "f4", ("y", "x"), **comp)
                var.units = units
                var.long_name = long_name
                if comment:
                    var.comment = comment
                var[:] = result[name]

save_timeseries_netcdf 

save_timeseries_netcdf(result: Dict, filepath: Union[str, Path], include_labels: bool = True, compression_level: int = 4) -> None

Save EQ-INSAR time series to NetCDF format.

Creates a NetCDF file with the time dimension for time series data, suitable for ML training datasets.

PARAMETER DESCRIPTION
result

Output from generate_timeseries

TYPE: dict

filepath

Output file path

TYPE: str or Path

include_labels

Include segmentation labels

TYPE: bool DEFAULT: True

compression_level

Compression level (0-9)

TYPE: int DEFAULT: 4

Source code in src/eq_insar/io/netcdf.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
def save_timeseries_netcdf(
    result: Dict,
    filepath: Union[str, Path],
    include_labels: bool = True,
    compression_level: int = 4,
) -> None:
    """
    Save EQ-INSAR time series to NetCDF format.

    Creates a NetCDF file with the time dimension for time series data,
    suitable for ML training datasets.

    Parameters
    ----------
    result : dict
        Output from generate_timeseries
    filepath : str or Path
        Output file path
    include_labels : bool
        Include segmentation labels
    compression_level : int
        Compression level (0-9)
    """
    if not _check_netcdf4():
        raise ImportError(
            "netCDF4 is required for NetCDF export. "
            "Install with: pip install netCDF4"
        )

    import netCDF4

    filepath = Path(filepath)
    filepath.parent.mkdir(parents=True, exist_ok=True)

    X_km = result["X_km"]
    Y_km = result["Y_km"]
    meta = result.get("metadata", {})
    timeseries = result["timeseries"]

    nt, ny, nx = timeseries.shape
    n_pre = meta.get("n_pre", 0)
    n_event = meta.get("n_event", 1)
    n_post = meta.get("n_post", 0)

    with netCDF4.Dataset(filepath, "w", format="NETCDF4") as nc:
        # Global attributes
        nc.title = "EQ-INSAR Synthetic InSAR Time Series"
        nc.institution = "EQ-INSAR Package"
        nc.source = "Synthetic earthquake deformation model"
        nc.history = f"Created {datetime.now().isoformat()}"
        nc.Conventions = "CF-1.8"

        # Time series metadata
        nc.n_pre_event = n_pre
        nc.n_event = n_event
        nc.n_post_event = n_post
        nc.output_type = meta.get("output_type", "phase")

        # Source parameters
        nc.earthquake_Mw = meta.get("Mw", np.nan)
        nc.earthquake_depth_km = meta.get("depth_km", np.nan)
        nc.earthquake_strike_deg = meta.get("strike_deg", np.nan)
        nc.earthquake_dip_deg = meta.get("dip_deg", np.nan)
        nc.earthquake_rake_deg = meta.get("rake_deg", np.nan)

        # Create dimensions
        nc.createDimension("time", nt)
        nc.createDimension("y", ny)
        nc.createDimension("x", nx)

        # Coordinates
        time_var = nc.createVariable("time", "i4", ("time",))
        time_var.units = "frame number"
        time_var.long_name = "Time frame index"
        time_var[:] = np.arange(nt)

        x_var = nc.createVariable("x", "f4", ("x",))
        x_var.units = "km"
        x_var.long_name = "Easting"
        x_var[:] = X_km[0, :]

        y_var = nc.createVariable("y", "f4", ("y",))
        y_var.units = "km"
        y_var.long_name = "Northing"
        y_var[:] = Y_km[:, 0]

        # Frame type indicator
        frame_type = nc.createVariable("frame_type", "i1", ("time",))
        frame_type.long_name = "Frame type (0=pre, 1=event, 2=post)"
        frame_type.flag_values = np.array([0, 1, 2], dtype=np.int8)
        frame_type.flag_meanings = "pre_event event post_event"
        types = np.zeros(nt, dtype=np.int8)
        types[n_pre:n_pre + n_event] = 1
        types[n_pre + n_event:] = 2
        frame_type[:] = types

        # Compression settings
        comp = {"zlib": True, "complevel": compression_level}

        # Time series data
        ts_var = nc.createVariable("timeseries", "f4", ("time", "y", "x"), **comp)
        ts_var.units = "rad" if meta.get("output_type") == "phase" else "m"
        ts_var.long_name = "InSAR time series"
        ts_var[:] = timeseries

        # Labels
        if include_labels and "labels" in result:
            label_var = nc.createVariable("labels", "i1", ("time", "y", "x"), **comp)
            label_var.long_name = "Deformation mask (1=deformation, 0=no deformation)"
            label_var[:] = result["labels"]

        # Static fields (not time-varying)
        if "los_displacement" in result:
            los_var = nc.createVariable("los_displacement", "f4", ("y", "x"), **comp)
            los_var.units = "m"
            los_var.long_name = "LOS displacement (static, no noise)"
            los_var[:] = result["los_displacement"]

load_netcdf 

load_netcdf(filepath: Union[str, Path]) -> Dict

Load EQ-INSAR data from NetCDF file.

PARAMETER DESCRIPTION
filepath

Path to NetCDF file

TYPE: str or Path

RETURNS DESCRIPTION
result

Dictionary with data arrays and metadata

TYPE: dict

Source code in src/eq_insar/io/netcdf.py 
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
def load_netcdf(filepath: Union[str, Path]) -> Dict:
    """
    Load EQ-INSAR data from NetCDF file.

    Parameters
    ----------
    filepath : str or Path
        Path to NetCDF file

    Returns
    -------
    result : dict
        Dictionary with data arrays and metadata
    """
    if not _check_netcdf4():
        raise ImportError(
            "netCDF4 is required for NetCDF import. "
            "Install with: pip install netCDF4"
        )

    import netCDF4

    filepath = Path(filepath)

    result = {}
    metadata = {}

    with netCDF4.Dataset(filepath, "r") as nc:
        # Load coordinates
        x_km = nc.variables["x"][:]
        y_km = nc.variables["y"][:]
        X_km, Y_km = np.meshgrid(x_km, y_km)
        result["X_km"] = X_km
        result["Y_km"] = Y_km

        # Load data variables
        for var_name in nc.variables:
            if var_name not in ["x", "y", "time", "frame_type"]:
                result[var_name] = nc.variables[var_name][:]

        # Load global attributes as metadata
        for attr in nc.ncattrs():
            value = nc.getncattr(attr)
            if attr.startswith("earthquake_"):
                key = attr.replace("earthquake_", "")
                metadata[key] = value
            else:
                metadata[attr] = value

    result["metadata"] = metadata
    return result