Good afternoon,
— Context
I am using otb within a self-made python class in order to call for various otb functions easily within a chain of operations. Among those i use the otb Haralick texture extraction which works fine. However the output seems a bit “off” (cf pictures)
Isolated pixels, but most importantly wide areas have an output value for the different haralick indices set to 0, which is the “no data value” of my file. The weirdest is that, the “patches” all have the same values across the all haralick indices. while the “isolated” pixels all behave differently.
I wonder if this is a normal behaviour of the otb haralick texture extraction.
Thank you for your help,
Jules Michard,
Regards.
—Parameters for OTB haralcik texture extraction
- xrand and yrad 1
- number of bins 16
- offset default
- haralick indices “simple”
- min and max are set to the min/max +/- 2*std of the image
I work with OTB 7.1 and can’t update it
— Code
################## IMPORTS
import os
import time
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import numpy as np
import rasterio
##### Faire l'import d'otb à la fin sinon tout plante
import otbApplication as otb
#------------------A modifier
fname="corona_30061966"
#--------------------------- CHEMINS --------------------------
class texture:
def __init__(self,filename):
year=filename[-4:]
self.path_corona="/home/jovyan/data_agata/"+str(year)+"/Images/"+str(filename)+"/"
self.path_mask="/home/jovyan/data_agata/"+str(year)+"/Images/"+str(filename)+"/" +str(filename)
self.path_mnt='/home/jovyan/data_agata/MNT'
self.path_indices = "/home/jovyan/data_agata/"+str(year)+"/Indices/"
self.raster_fname=str(filename)
self.corrected_raster=self.path_corona+str(filename)+'.tif'
#------------------ Calculs des différents indices sur l'image CORONA (composée et clippé) ------------------
def bin_val(self,fenetre):
if fenetre<=10:
return int(64)
else:
return int(16)
#Bande: renseigner 1, 2 ,3 ou 4 pour les bandes 4,5,6, ou è respectivement
#fenetre: int
def haralick(self,filename,bande,nom_bande,fenetre):
entree= self.corrected_raster
#sortie A changer pour aboutir dans dossier INDICES + ajouter taille radius
sortie=self.path_indices + str(filename)+'_haralick_{}.tif'.format(fenetre)
app=None
app = otb.Registry.CreateApplication("HaralickTextureExtraction")
#Entrée
app.SetParameterString("in", entree)
#Bande sur laquelle c'est calculé
app.SetParameterInt("channel",bande)
#Borner les indices d'Haralick (extraire max er min)
bande_objet=rasterio.open(entree).read(bande)
app.SetParameterInt("parameters.nbbin",self.bin_val(fenetre=fenetre))
app.SetParameterFloat("parameters.min",float(bande_objet.mean()-2*bande_objet.std()))
app.SetParameterFloat("parameters.max",float(bande_objet.mean()+2*bande_objet.std()))
#Taille de la fenêtre
app.SetParameterInt("parameters.xrad", fenetre)
app.SetParameterInt("parameters.yrad", fenetre)
#Offset, c'est quoi ?
#Choix des indices ---> Lequel choisir
app.SetParameterString("texture","simple")
app.SetParameterString("out", sortie)
app.ExecuteAndWriteOutput()
#On réaligne les étendues si décalage
ds_ref=gdal.Open(entree)
gt_ref=ds_ref.GetGeoTransform
ulx, xres, xskew, uly, yskew, yres = ds_ref.GetGeoTransform()
lrx = ulx + (ds_ref.RasterXSize * xres)
lry = uly + (ds_ref.RasterYSize * yres)
#unset nodata ne marche pas
commande="/opt/anaconda/envs/env_otb/bin/gdal_edit.py -a_srs epsg:32718 -a_ullr {ulx} {uly} {lrx} {lry} {sortie}".format(ulx=ulx,uly=uly,lrx=lrx,lry=lry,sortie=str(sortie))
os.system(commande)
ds_ref=None
return
###############Execution
def compute_index(self,filename):
win_sizes=[1]
dataset = rasterio.open(self.corrected_raster)
#On a 1 bande sur CORONA
for n_band in range(1,dataset.count+1):
print('Pour la bande numéro {n_band}/{tot}'.format(n_band=n_band,tot=dataset.count))
for fen in win_sizes:
print("haralick pour",fen)
self.haralick(filename,n_band,'',fen) #2 min avec hlick_rad = 11
calcul=texture(fname)
calcul.compute_index(fname)