keeping the coordinate system of raster files in the resulting raster file after operation with numpy

Question

I'm trying to read as array two raster files using ReadAsArray in GDAL/python, perform some map algebra operation on it using numpy and then write the resulting raster file in GTiff format using WriteArray in GDAL/Python. The spatial reference system of the input raster files are the same. I want to know the resulting GTiff file will have the same coordinate system or not?

Answer 1

The only Python WriteArray() method I know of is gdal.Band.WriteArray() which is called on an existing file. It would be at the point where the file is created that it's coordinate system and geotransform is established. So a simplistic use of gdal.Driver.Create() and gdal.Band.WriteArray() will not preserve the coordinate system.

You have a couple options.

1. Create the output file with gdal.Driver.Create() and carefully set the projection and geotransform to match the source image.
2. Alternatively write the array using the gdal_array.SaveArray() method which includes a prototype input file from which the georeferencing info should be closed.

eg. (not actually tested)

img1 = gdal_array.LoadFile("abc.tif")
img2 = img1 * 2
gdal_array.SaveArray(img2, "out.tif", "GTiff", gdal.Open("abc.tif"))