everyone, I'm looking to integrate a bash command into my python code to calculate indices. My problem is that I want to have an output image with a band for each of the calculated indices, but I can't integrate these indices by the bash command into my 'im_index' matrix created with my python code. I don't see how to link both of them... Do you have any idea?
import numpy as np
import sys
import os
import spectral as sp
from scipy import ndimage
import pylab as pl
from math import *
import spectral.io.envi as envi
#------------------------------------
def reject_outliers(data, m=1):
return data[abs(data - np.mean(data)) < m * np.std(data)]
#------------------------------------
def find_nearest(array, value):
#For a given value, find the nearest value in an array
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
#------------------------------------
#Open existing dataset
src_directory = "/d/afavro/Bureau/4_reflectance/"
dossier = os.listdir (src_directory)
print(dossier)
for fichier in dossier:
print (fichier)
ssrc_directory = "/d/afavro/Bureau/4_reflectance/" + fichier + "/"
rasters = os.listdir (ssrc_directory)
print(rasters)
OUTPUT_FOLDER = "/d/afavro/Bureau/5_indices2/" + 'indices_' + fichier + '/'
print(OUTPUT_FOLDER)
if not os.path.exists(OUTPUT_FOLDER):
os.makedirs(OUTPUT_FOLDER)
for image in rasters:
print(image)
name, ext = os.path.splitext(image)
if ext == '.hdr':
img = sp.open_image(ssrc_directory + image)
print(image)
im_HS = img[:,:,:]
cols = im_HS.shape[0] # Number of column
rows = im_HS.shape[1] # Number of lines
bands = im_HS.shape[2] # Number of bands
NbPix = cols * rows # Number of pixels
#Get wavelengths from hdr file
wv = np.asarray(img.bands.centers)
if len(wv) == 0 :
print("Wavelengths not defined in the hdr file")
sys.exit("Try again!")
if wv[0] > 100:
wv=wv*0.001 # Convert to micrometers if necessary
im_HS=im_HS.reshape(NbPix, bands)
#Compute HC index------------------------------------------------------
Nind=4 # Number of indice to be computed
im_index=np.zeros((cols*rows, Nind))
names = []
##NDVI computation
names.append('NDVI')
bande_ref=[0.67, 0.8]
bRef0 = find_nearest(wv,bande_ref[0])
bRef1 = find_nearest(wv,bande_ref[1])
#Check if the required specral bands are available
if (np.abs(wv[bRef0]-bande_ref[0])<=0.1 and np.abs(wv[bRef1]-bande_ref[1])<=0.1):
b0 = im_HS[:, bRef0]
b1 = im_HS[:, bRef1]
index = (b0 - b1) / (b0 + b1)
else:
index = np.zeros_like(im_HS[:,0])
print("Wavelengths selection problem, NDVI not computed")
im_index[:,0]= index
# bash command :
inRaster = ssrc_directory + image
print(inRaster)
outRaster = OUTPUT_FOLDER + 'indices_' + image
print (outRaster)
cmd = 'otbcli_RadiometricIndices -in inRaster -list Soil:BI Vegetation:MSAVI Vegetation:SAVI -out outRaster'
os.system(cmd)
#saving
im_index=im_index.reshape(cols, rows, Nind)
file_image = OUTPUT_FOLDER + "indices2_" + fichier
header = envi.read_envi_header(ssrc_directory + image)
header ['description'] = "fichier d'origine " + image
header ['band names'] = ['NDVI', 'Sober filter', 'NDWI', 'IB(1)', 'IB(2)']
del header['wavelength units']
del header['wavelength']
sp.envi.save_image(file_image + '.hdr', im_index, metadata=header, force = True, interleave = 'bsq')
Assuming this is the code you are actually asking about:
inRaster = ssrc_directory + image
print(inRaster)
outRaster = OUTPUT_FOLDER + 'indices_' + image
print (outRaster)
cmd = 'otbcli_RadiometricIndices -in inRaster -list Soil:BI Vegetation:MSAVI Vegetation:SAVI -out outRaster'
os.system(cmd)
Of course, inRaster inside of singe quotes is just a literal string; to interpolate the variable's value you can say
cmd = 'otbcli_RadiometricIndices -in ' + inRaster + \
' -list Soil:BI Vegetation:MSAVI Vegetation:SAVI -out ' + \
outRaster
or
cmd = 'otbcli_RadiometricIndices -in {0} -list Soil:BI Vegetation:MSAVI Vegetation:SAVI -out {1}'.format(
inRaster, outRaster)
or a number of other string interpolation techniques in Python (legacy % formatting, f-string, etc). But a better solution is to replace os.system with the more flexible and versatile subprocess, as suggested even in the os.system documentation.
subprocess.run([
'otbcli_RadiometricIndices',
'-in', inRaster,
'-list', 'Soil:BI', 'Vegetation:MSAVI', 'Vegetation:SAVI',
'-out', outRaster], check=True)
subprocess.run was introduced in Python 3.5; if you need compatibility with older versions, try subprocess.check_call or even the crude subprocess.call.