ENH : Python script to compute average spectra in range and in azimuth for a Tiff image

python_src/utils/compute_spectra_ImageTiff.py

+# coding:utf-8
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import argparse
+try:
+    import gdal
+except ImportError:
+    import osgeo.gdal as gdal
+def calcul_spectres_tiff(ficima,output_path, output_filename_az, output_filename_dist):
+    """
+    Compute the average spectra in range and in azimuth from a complex SAR image (mandatory format .tif)
+    """
+    # Average spectrum in azimut
+    ###################################
+
+    print("Compute the average spectrum in azimut after loading the modules python and gdal")
+    #Read the image tiff
+    zone = read_tiff(ficima)
+    print ("Dimensions of the tif file ",np.shape(zone))
+    (nb_lig, nb_col) = np.shape(zone[:,:,0])
+    print ("nb lig = ",nb_lig," nb col = ", nb_col)
+
+    lig = np.zeros(nb_lig, dtype = complex)
+    spectre_az = np.zeros(nb_lig, dtype = 'float')
+
+    for i in range(nb_col):
+        lig = zone[:,i,0] + 1j*zone[:,i,1] #All the lines in complex
+        spectre_az = spectre_az + np.abs(np.fft.fft(lig))**2
+
+    spectre_az = np.sqrt(spectre_az) / float(nb_col)
+
+    #Display and Save the spectrum
+    display_spectrum(spectre_az, "Average spectrum in azimut", output_path, output_filename_az)
+
+    # Compute the average spectrum in range
+    #####################################
+
+    print("Compute the average spectrum in range")
+    col = np.zeros(nb_col, dtype = complex)
+    spectre_dist = np.zeros(nb_col, dtype = 'float')
+
+    for i in range(nb_lig):
+        col = zone[i,:,0] + 1j*zone[i,:,1] #All the columns in complex
+        spectre_dist = spectre_dist + np.abs(np.fft.fft(col))**2
+
+    spectre_dist = np.sqrt(spectre_dist) / float(nb_lig)
+
+    #Display and Save the spectrum
+    print("Display and Save the spectrum")
+    display_spectrum(spectre_dist, "Average spectrum in range", output_path, output_filename_dist)
+
+    plt.show()
+
+def read_tiff(filename):
+    """
+    Create a numpy array from the input image tiff.
+    """
+    ds = gdal.Open(filename)
+
+    nb_lig = ds.RasterYSize
+    nb_col = ds.RasterXSize
+    nb_bands = ds.RasterCount
+
+    zone = np.zeros((nb_lig, nb_col, nb_bands))
+
+    for i in range(nb_bands):
+        zone[:,:,i] = np.array(ds.GetRasterBand(i+1).ReadAsArray())
+
+    return zone
+
+def display_spectrum(spectre, title, path, filename):
+    """
+    Display the normalized spectrum.
+    """
+    nb_val = np.size(spectre)
+    xvals = np.arange(nb_val) / float(nb_val)
+
+    plt.figure()
+    plt.plot(xvals, spectre)
+    plt.ylim(0, 1.1*max(spectre))
+    plt.title(title + "\n" + filename)
+    plt.savefig(os.path.join(path,filename))
+    plt.draw()
+
+
+if __name__ == "__main__":
+    description = "Compute of the spectra of the image tif"
+    parser = argparse.ArgumentParser(description=description)
+    parser.add_argument("path_input", type=str, help="Path of the input Tiff Image to get its spectrum. There is not its extension in the path")
+    parser.add_argument("path_output_dir", type=str, help="Path of the output directory where the spectra of the input will be saved")
+    args = parser.parse_args()
+
+
+    print("Usage : python compute_spectra_ImageTiff.py <path_input_tif_image without the extension file> <path_output_directory>")
+    print("The Tiff Image must have 2 bands in float/int. The script cannot read the image with one band in complex")
+    path_input = args.path_input
+    file_input = path_input.split('/')[-1]
+    input_path = os.path.dirname(path_input)
+    if(input_path==''):
+        input_path = '.'
+    nomima = file_input.split('.')[0]
+    ficima=os.path.join(input_path,nomima+".tif")
+    output_path=args.path_output_dir
+    output_filename_az="spectrum_az_"+nomima
+    output_filename_dist="spectrum_dist_"+nomima
+
+    calcul_spectres_tiff(ficima,output_path, output_filename_az, output_filename_dist)