Process Satellite Subsidence Data with QGIS and Python on Linux (2025)

Geospatial developers working with Interferometric Synthetic Aperture Radar (InSAR) data from missions like ISRO's or NASA satellites face a common challenge: converting raw subsidence measurements into actionable visualizations and analysis. Recent US-Indian space missions mapping extreme subsidence in Mexico City have generated massive SAR datasets that require specialized processing pipelines.

This guide walks you through setting up a complete workflow to process satellite subsidence data on Linux using QGIS as your visualization layer and Python for data manipulation.

Why QGIS + Python for Subsidence Analysis

When you're working with InSAR data from space missions, you need tools that handle:

• Georeferenced raster datasets (satellite imagery)
• Large NetCDF or HDF5 elevation change files
• Custom calculations on grid data
• Publication-quality map outputs

QGIS provides the GUI interface for quick exploration, while Python automates repetitive processing tasks and handles batch operations across multiple datasets.

Prerequisites and Setup

First, ensure your Linux system has the necessary geospatial libraries installed:

sudo apt-get update
sudo apt-get install qgis python3-pip gdal-bin proj-bin
pip3 install rasterio fiona numpy scipy matplotlib geopandas xarray netCDF4

For Ubuntu 22.04 LTS and later, this stack runs reliably with Python 3.10+. Verify your GDAL installation supports the satellite data formats you're working with:

gdalinfo --formats | grep -E "HDF5|NetCDF"

Step 1: Download and Organize Raw SAR Data

Satellite missions typically provide subsidence data in one of two formats:

• GeoTIFF files (easiest to work with)
• NetCDF or HDF5 hierarchical formats (common from research institutions)

Create a project directory structure:

mkdir -p ~/subsidence_project/{raw_data,processed,output}
cd ~/subsidence_project

Once downloaded, inspect the data metadata:

import rasterio
import xarray as xr

# For GeoTIFF files
with rasterio.open('raw_data/subsidence_mexico_city.tif') as src:
    print(f"CRS: {src.crs}")
    print(f"Bounds: {src.bounds}")
    print(f"Shape: {src.height} x {src.width}")
    print(f"Data type: {src.dtypes[0]}")

# For NetCDF files
ds = xr.open_dataset('raw_data/sar_displacement.nc')
print(ds)

Step 2: Pre-process Subsidence Data with Python

Raw InSAR data often requires unit conversion, projection adjustments, and artifact removal. This Python script handles common preprocessing tasks:

import rasterio
from rasterio.transform import from_bounds
from rasterio.crs import CRS
import numpy as np
from scipy import ndimage
import warnings
warnings.filterwarnings('ignore')

def preprocess_subsidence(input_path, output_path, units='mm'):
    """
    Preprocess raw subsidence data:
    - Convert units (mm/cm/meters)
    - Apply spatial filter to remove noise
    - Generate quality mask
    """
    with rasterio.open(input_path) as src:
        data = src.read(1).astype(float)
        profile = src.profile
        
        # Handle nodata values
        mask = data == src.nodata
        data[mask] = np.nan
        
        # Convert to millimeters if needed
        if units == 'cm':
            data *= 10
        elif units == 'meters':
            data *= 1000
        
        # Apply Gaussian filter to reduce speckle noise
        # Preserve edges with bilateral-like approach
        filtered = ndimage.gaussian_filter(data, sigma=1.5)
        
        # Update profile for output
        profile.update(dtype=rasterio.float32, nodata=np.nan)
        
        # Write processed data
        with rasterio.open(output_path, 'w', **profile) as dst:
            dst.write(filtered.astype(rasterio.float32), 1)
    
    return output_path

# Process your subsidence data
preprocess_subsidence(
    'raw_data/subsidence_mexico_city.tif',
    'processed/subsidence_clean.tif',
    units='cm'
)

Step 3: Load Processed Data into QGIS

Now open QGIS and load your processed raster:

1. Launch QGIS: qgis processed/subsidence_clean.tif
2. Add basemap for context:
   • Go to Browser → XYZ Tiles → Add OpenStreetMap
3. Style the subsidence layer:
   • Right-click layer → Properties → Symbology
   • Choose Singleband Pseudocolor renderer
   • Select Turbo or Spectral colormap
   • Set Min/Max values based on your subsidence range (e.g., -500mm to +100mm)

Step 4: Extract Statistics by Administrative Boundary

To analyze subsidence by region (like Mexico City boroughs), use vector data and zonal statistics:

import geopandas as gpd
from rasterio.mask import mask
import rasterio

# Load administrative boundaries (download from INEGI)
admin_bounds = gpd.read_file('raw_data/mexico_city_boroughs.geojson')

# Open subsidence raster
with rasterio.open('processed/subsidence_clean.tif') as src:
    stats = []
    
    for idx, row in admin_bounds.iterrows():
        geom = [row.geometry]
        borough_name = row['nombre']
        
        try:
            # Clip raster to borough
            clipped, transform = mask(
                src, geom, crop=True, nodata=np.nan
            )
            values = clipped[0][~np.isnan(clipped[0])]
            
            if len(values) > 0:
                stats.append({
                    'borough': borough_name,
                    'mean_subsidence_mm': np.mean(values),
                    'max_subsidence_mm': np.max(values),
                    'pixel_count': len(values)
                })
        except Exception as e:
            print(f"Error processing {borough_name}: {e}")
    
    results_df = pd.DataFrame(stats)
    results_df.to_csv('output/subsidence_by_borough.csv', index=False)
    print(results_df)

Comparison: QGIS vs ArcGIS for Subsidence Analysis

| Feature | QGIS (Linux) | ArcGIS Pro (Windows/Mac) | |---------|--------------|-------------------------| | Cost | Free/Open Source | $2,000/year | | Python Integration | PyQGIS (built-in) | ArcPy (licensing required) | | Raster Processing Speed | Good for <500MB files | Excellent for large datasets | | Linux Support | Native | Not officially supported | | InSAR Tool Libraries | Community QGIS plugins | Limited commercial plugins | | Learning Curve | Moderate | Steep |

For developers processing subsidence data on Linux servers, QGIS + Python is the clear winner due to native support and no licensing barriers.

Common Pitfalls and Solutions

Projection Mismatches: Always verify CRS matches your administrative boundaries. Mexico City data should use EPSG:32614 (UTM Zone 14N) or reproject to EPSG:4326 (WGS84).

Memory Issues with Large Rasters: Process in tiles using rasterio.windows.Window() instead of loading entire arrays.

Nodata Handling: Ensure your filtered output preserves original nodata values; don't convert them to zeros, as this skews statistics.

Next Steps

Once you've processed subsidence data locally:

1. Publish to web: Load results into Geoserver/QGIS Server
2. Time-series analysis: Stack multiple mission datasets and track subsidence over years
3. Machine learning: Train models to predict future subsidence based on groundwater extraction patterns

With this workflow, you can process satellite subsidence data independently without proprietary GIS software, perfect for researchers and developers working with freely available mission data.