Dear community, i'm new to gdal and opencv and I would like to run bilateral filter operation on a geotiff with up to 77 channels. So far i managed to work out a version handling a 1 layered image and performing the filter operation. Geoinformation is handed over by saving the image two times. first time opencv filter result is written to image that is opened again with gdal and saved with geoinformation from former source image. Is there away to directly hand over the image data to gdal without this two times export routine? How can I handle multiple layers? How can I make GDAL imports operational for opencv filter operations? Greetings confugere

import cv2 import numpy as np from cv2 import waitKey from osgeo import gdal, osr

read geotiff

def readFile(filename): filehandle = gdal.Open(filename) band1 = filehandle.GetRasterBand(1) geotransform = filehandle.GetGeoTransform() geoproj = filehandle.GetProjection() band1data = band1.ReadAsArray() xsize = filehandle.RasterXSize ysize = filehandle.RasterYSize return xsize,ysize,geotransform,geoproj,band1data

write geotiff

def writeFile(filename,geotransform,geoprojection,data): xy = data.shape format = "GTiff" driver = gdal.GetDriverByName(format) dst_datatype = gdal.GDT_Byte dst_ds = driver.Create(filename,xy[1],xy[0],1,dst_datatype) dst_ds.GetRasterBand(1).WriteArray(data) dst_ds.SetGeoTransform(geotransform) dst_ds.SetProjection(geoprojection) return 1

def array_to_raster(array,x_pixels,y_pixels,xmin,ymax,psize,proje): """Array > Raster Save a raster from a C order array.

:param array: ndarray
"""
dst_filename = 'D:/testgdal.tif'


# You need to get those values like you did.
x_pixels = x_pixels  # number of pixels in x
y_pixels = y_pixels  # number of pixels in y
PIXEL_SIZE = psize  # size of the pixel...        
x_min = xmin  
y_max = ymax  # x_min & y_max are like the "top left" corner.
wkt_projection = proje

driver = gdal.GetDriverByName('GTiff')

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

dataset.SetGeoTransform((
    x_min,    # 0
    PIXEL_SIZE,  # 1
    0,                      # 2
    y_max,    # 3
    0,                      # 4
    -PIXEL_SIZE))  

dataset.SetProjection(wkt_projection)
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.

#

#

test bilateral filter

datafile = gdal.Open("D:/satellitendaten/test.tif") img = cv2.imread("D:/satellitendaten/test.tif")

bilateralResult = cv2.bilateralFilter(img,9,50,50)

filter does not work with thisone:

imgdal = cv2.imread("D:/satellitendaten/test.tif", cv2.IMREAD_ANYDEPTH)

#

#

blurnum = np.asarray(blur)

cv2.imshow("salida",blur)

cv2.imwrite("D:/satellitendaten/out.tif",bilateralResult)

#

#

read geoinformation from source and save it to new output

source

(xtest,ytest,transformtest,projtest,datatest) = readFile("D:/satellitendaten/test.tif")

output

(xblur,yblur,transformblur,projblur,bilateraldata) = readFile("D:/satellitendaten/out.tif")

print ("Ergebnis test: ",xtest," ",ytest," ", transformtest, " ", projtest, " ") print ("Ergebnis blur: ",xblur," ",yblur," ", transformblur, " ", projblur, " ") print transformtest[3]

writeFile("D:/satellitendaten/outgeo.tif",transformtest,projtest,bilateraldata)