GeotiffOffsets_2bands.py

#ISCE Offset outputs to Geotiff
#Written By Saeid Aminjafari & Mohammadnia
#Email: saeed.aminjafari@natgeo.su.se
#
#This code extracts the bands of a multi-band image generated by ISCE software (such as Range and Azimuth offsets) and saves the result as a Geotiff image.
#To run the code in terminal, change the desired fields based on comments and type: >>python3 GeotiffOffsets_2bands.py
#Your result will be stored in current directory.

from shutil import copyfile
from osgeo import gdal            ## GDAL support for reading virtual files
import os                         ## To create and remove directories
import matplotlib.pyplot as plt   ## For plotting
import numpy as np                ## Matrix calculations
import glob
import gdal
from gdalconst import *
from osgeo import osr

def array_to_raster(array):

    dst_filename = InputName + '_rg.tiff'        #Range offset output name
    #dst_filename = InputName + '_amp.tiff'      #Azimuth offset output name

    driver = gdal.GetDriverByName('GTiff')

    dataset = driver.Create(
        dst_filename,
        x_pixels,
        y_pixels,
        1,
        gdal.GDT_Float32, )

    dataset.SetGeoTransform((
        x_min,
        PIXEL_SIZE,
        0,
        y_max,
        0,
        -PIXEL_SIZE))  

    dataset.SetProjection(prj_wkt)
    dataset.GetRasterBand(1).WriteArray(array)
    dataset.FlushCache()  # Write to disk.
    return dataset, dataset.GetRasterBand(1)  #If you need to return, remenber to return  also the dataset because the band don`t live without dataset. 

InputName = 'filt_topophase.unw.geo'                           #Change the file name
ds = gdal.Open(InputName, gdal.GA_ReadOnly)                    #Open your file
x_pixels = ds.RasterXSize                                      #number of pixels in x direction
y_pixels = ds.RasterYSize                                      #number of pixels in y direction
PIXEL_SIZE = ds.GetGeoTransform()[1]                           #Pixel size
x_min = ds.GetGeoTransform()[0]                                #Upper left corner (X)
y_max = ds.GetGeoTransform()[3]                                #Upper left corner (y)
prj_wkt = ds.GetProjectionRef()                                #Projection
rg = ds.GetRasterBand(1).ReadAsArray()                         #Range part of the data
az = ds.GetRasterBand(2).ReadAsArray()                         #Azimuth part of the data
array_to_raster(rg)                                            #Saving range offsets to Geotiff
#array_to_raster(az)                                            #Saving azimuth offsets to Geotiff