How to save raster with projection into Python

Question

I applied a raster processing that is already in a known coordinate system with this code,

import cv2
import numpy as np
from matplotlib import pyplot as src
img = cv2.imread(r'E:2円_PROJETS_DISK_E\threshold621円.tif',0)
# Otsu's thresholding after Gaussian filtering
blur = cv2.GaussianBlur(img,(5,5),0)
ret1,th1 = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
for i in xrange(1):
 plt.imshow(images[i*3+2],'gray')
 plt.xticks([]), plt.yticks([])
plt.show()
cv2.imwrite( r'E:2円_PROJETS_DISK_E\threshold621円-1.tif',th1);

but when I save my raster at the end of this script, I get a non-georeferenced TIFF raster

How can I keep the same coordinate system of the initial raster (without tranforming the output raster into local coordinates)

Since I am new to Python, and I have no knowledge of Python, i would like someone to correct me my script please.

• gene: You cannot georeference a raster with Opencv. You need to use a geospatial module as GDAL or rasterio
• hamza: i have installed GDAL but i dont know how to use the command, can you give me the correcte script please?
• gene: Transforming a image to a shape defined by image corners in earth coordinates - Python for example

Answer 1

thank you for the answer, but it doesnt works for me

i tried this and it works, hope it will help another:

import cv2
import numpy as np
import gdal
in_imgpath = r'E:2円_PROJETS_DISK_E\Raster2.tif'
img = cv2.imread(in_imgpath ,0)
dataset1 = gdal.Open(in_imgpath)
projection = dataset1.GetProjection()
geotransform = dataset1.GetGeoTransform()
# Otsu's thresholding after Gaussian filtering
blur = cv2.GaussianBlur(img,(5,5),0)
ret1,th1 = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
kernal = np.ones((3,3), np.uint8)
dilation = cv2.dilate(th1, kernal, iterations=2)
erosion = cv2.erode(dilation, kernal, iterations=1)
opening = cv2.morphologyEx(erosion, cv2.MORPH_OPEN, kernal, iterations=3)
closing = cv2.morphologyEx(opening, cv2.MORPH_CLOSE, kernal, iterations=4)
out_imgpath = r'E:2円_PROJETS_DISK_E\Raster2.tif'
cv2.imwrite(out_imgpath ,closing)
dataset2 = gdal.Open(out_imgpath, gdal.GA_Update)
dataset2.SetGeoTransform( geotransform )
dataset2.SetProjection( projection )

Answer 2

You should be able to set the projection and the geotransform of your new raster using the information of the previous raster.

Can you test the following code? (You can append it to your code or create a new script a run it after the code you provided).

import gdal
# open original raster and get projection and geotransform
original_ds = gdal.Open(r'E:2円_PROJETS_DISK_E\threshold621円.tif', 0)
sr = ds.GetProjection()
gt = ds.GetGeoTransform()
del original_ds
# open target raster (in editing mode) and set the projection and geotransform
ds = gdal.Open(r'E:2円_PROJETS_DISK_E\threshold621円-1.tif', 1)
ds.SetProjection(sr)
ds.SetGeoTransform(gt)
# save and close target raster
del ds

Answer 3

You can read the image using gdal.ReadAsArray() and write the image with gdal as well.

import cv2
import numpy as np
import gdal
in_imgpath = r'E:2円_PROJETS_DISK_E\Raster2.tif'
dataset1 = gdal.Open(in_imgpath, gdal.GA_ReadOnly)
projection = dataset1.GetProjection()
geotransform = dataset1.GetGeoTransform()
img = dataset1.ReadAsArray()
rows, cols = img.shape # assuming it's a 2D array
dataset1 = None
# Otsu's thresholding after Gaussian filtering
blur = cv2.GaussianBlur(img,(5,5),0)
ret1,th1 = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
kernal = np.ones((3,3), np.uint8)
dilation = cv2.dilate(th1, kernal, iterations=2)
erosion = cv2.erode(dilation, kernal, iterations=1)
opening = cv2.morphologyEx(erosion, cv2.MORPH_OPEN, kernal, iterations=3)
closing = cv2.morphologyEx(opening, cv2.MORPH_CLOSE, kernal, iterations=4)
out_imgpath = r'E:2円_PROJETS_DISK_E\Raster2.tif'
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(out_imgpath, cols, rows, n, gdal.GDT_Byte)
dataset.SetGeoTransform(geo_transform)
dataset.SetProjection(projection)
band.WriteArray(closing)
dataset = None
band = None

Answer 4

Just write the output array with GDAL along with the transform and projection information.

import cv2 
from osgeo import gdal
def writeArr2Tif(outFileName, arr_out, rows, cols, bandIndex, noData, geoTrans, proj):
 driver = gdal.GetDriverByName("GTiff")
 outDs = driver.Create(outFileName, rows, cols, bandIndex, gdal.GDT_Byte)
 outDs.SetGeoTransform(geoTrans)##sets same geotransform as input
 outDs.SetProjection(proj)##sets same projection as input
 outDs.GetRasterBand(bandIndex).WriteArray(arr_out)
 outDs.GetRasterBand(bandIndex).SetNoDataValue(noData)##if you want these values transparent
 outDs.FlushCache() ##saves to disk!!
 outDs = None
 print("done")
ds = gdal.Open(infile)
img_arr = ds.ReadAsArray()
width = ds.RasterXSize
height = ds.RasterYSize
nBand = ds.RasterCount 
trans = ds.GetGeoTransform()
proj = ds.GetProjection()
# Taking a matrix of size 5 as the kernel 
kernel = np.ones((3, 3), np.uint8) 
img_dilation = cv2.dilate(img_arr, kernel, iterations=1) 
writeArr2Tif(out_path_ero, img_erosion, height, width, 1, 255, trans, proj)