How to extract the phase band and convert to .grd or .tif format in order to plot interferogram by GMT/PyGMT or open in QGIS?

[leeyhgeo]

Hello,
I'm recently using ISCE (topsApp.py) to produce sentinel-1 interferograms and I would like to plot the outputs by GMT or PyGMT.
As a result, I converted the output file "filt_topophase.flat.geo" to .grd format by using gdal_translate -of GMT filt_topophase.flat.geo filt_topophase.flat.grd
However, I found that the converted .grd file seems to combine amplitude and phase bands together, rather than amplitude and phase separately.
Does anyone know how to extract the phase band only in order to plot the interferogram by GMT or PyGMT?

[riyasmjgeo]

Refer -b flag in gdal_translate : https://gdal.org/programs/gdal_translate.html
You can assign which band to export. I guess it could help.

[leeyhgeo]

Hello, thanks a lot for your advice. However, I checked the raster data and found that there is only one band. As a result, I cannot use gdal_translate -b to extract phase information from the interferogram.
I think it's a complex32 image according to the metadata below. Is it possible to separate the real and imaginary part of the image?

<VRTDataset rasterXSize="4653" rasterYSize="6066">
    <SRS>EPSG:4326</SRS>
    <GeoTransform>120.58458333333334, 0.00027777777777777696, 0.0, 24.97486111111111, 0.0, -0.0002777777777777779</GeoTransform>
    <VRTRasterBand dataType="CFloat32" band="1" subClass="VRTRawRasterBand">
        <SourceFilename relativeToVRT="1">filt_topophase.flat.geo</SourceFilename>
        <ByteOrder>LSB</ByteOrder>
        <ImageOffset>0</ImageOffset>
        <PixelOffset>8</PixelOffset>
        <LineOffset>37224</LineOffset>
    </VRTRasterBand>
</VRTDataset>

[riyasmjgeo]

Okay. I understood now. I don't know how to directly extract the complex part. However, here is an inefficient solution. Hopefully, this will work, but pls make sure the result is correct before proceeding.

#!/usr/bin/env python
from osgeo import gdal
import numpy as np

# Read the file and extract attributes 
input_path = './filt_topophase.flat.geo'
ds = gdal.Open(input_path, gdal.GA_ReadOnly)
transform = ds.GetGeoTransform()
proj = ds.GetProjection()

# Convert to numpy and extract the complex (interferogram) part
combined_prod = ds.GetRasterBand(1).ReadAsArray()
ifg = np.angle(combined_prod)
# Export the interferogram in GTiff
driver = gdal.GetDriverByName("GTiff")
driver.Register()
output = driver.Create("output.tif", xsize = ifg.shape[1], ysize = ifg.shape[0], bands=1, eType = gdal.GDT_Float64)
output.SetGeoTransform(transform)
output.SetProjection(proj)
outband = output.GetRasterBand(1)
outband.WriteArray(ifg)
outband.SetNoDataValue(np.nan)
outband.FlushCache()
outband = output = None

[leeyhgeo]

I just tried your code and it worked!!!
Thank you so much for this solution.

[ditafaith]

Hello @riyasmjgeo I just completed get the result after a pair processing.
However, I am confused about which file I should import to QGIS to see the displacement, such unwrapped and wrapped files.