1

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\threshold\621.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\threshold\621-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.

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 )
1

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\threshold\621.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\threshold\621-1.tif', 1)
ds.SetProjection(sr)
ds.SetGeoTransform(gt)

# save and close target raster
del ds

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.