Dear community I'd like to run bilateral filter operations on some Geotiffimages that can contain up to 70 layers. I'm new to opencv and gdal in python and so far managed to apply the bilateral filter on a test image with one layer and export it again as a geotiff. The export so far works by writing the image with opencv, open it with gdal again and write it again with gdal taking the geoinformation from the former input. I'd like to know how to handle multiple channels and to directly hand over the data imported with opencv to gdal and write it with the geoinformation gained from the source. Here is my code sofar:

''' Created on 24.11.2014

@author: Andreas ''' 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)