Geotiff TIFF Channels – How to Get Channels of a .tiff UAV Image and Maintain the Georeference

gdalgeotiff-tiffpythonrasterunmanned aerial vehicle

So, the idea is that I have a .tif UAV image. I need to get R,G,B,A channels into 4 separate .tif files by maintaining the georeferencing. So, for instance: Red.tif, Green.tif, Blue.tif and Alpha.tif. I am using this approach:

(red, green, blue, alpha) = np.transpose(img, axes = (2,0,1))

Next, I want to do some calculations with the channels. For instance:

result = ((red**2)+(blue**2))/(blue)

I use this code in order to make the result.tif…

import rasterio
with rasterio.open('path/to/Red_channel.tif') as f:
    red = f.read()
    profile = f.profile

with rasterio.open('path/to/Blue_channel.tif') as f:
    blue = f.read()

result = ((red**2)+(blue**2))/(blue) 

with rasterio.open('path/to/Output.tif', 'w', **profile) as dst:
    dst.write(result)

(source: Calculations with .tif images using matplotlib or rasterio)

Now, the problem is that when I try to use the Output.tif with gdalinfo to see the real coordinates, it does not show real coordinates, it shows pixel coordinates for each corner!! Any idea what is wrong and how I fix this?

Update:
I store each of the bands after:

(red, green, blue, alpha) = np.transpose(f, axes = (2,0,1))

like this:

with rasterio.open('path/to/Green.tif', 'w', **profile) as dst:
    dst.write(green)

and I have set:

profile["count"] = 4

This is the error I got:

Traceback (most recent call last):
  File "code.py", line 126, in bands
    (red, green, blue, alpha) = np.transpose(f, axes = (2,0,1))
  File "<__array_function__ internals>", line 5, in transpose
  File "/home/UbuntuUser/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 653, in transpose
    return _wrapfunc(a, 'transpose', axes)
  File "/home/UbuntuUser/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 55, in _wrapfunc
    return _wrapit(obj, method, *args, **kwds)
  File "/home/UbuntuUser/.local/lib/python3.8/site-packages/numpy/core/fromnumeric.py", line 44, in _wrapit
    result = getattr(asarray(obj), method)(*args, **kwds)
ValueError: axes don't match array

Best Answer

I tried the following with a GTiff image I split into R, G and B bands in QGIS:

wd = "C:\\Users\\Manuel\\Desktop\\test"
import rasterio, os
R = os.path.join(wd, "xR.tif")
B = os.path.join(wd, "xB.tif")
output = os.path.join(wd, "out.tif")

with rasterio.open(R) as f:
    red = f.read()
    profile = f.profile

with rasterio.open(B) as f:
    blue = f.read()

profile["dtype"] = "float64"
result = ((red**2)+(blue**2))/(blue) 

with rasterio.open(output, 'w', **profile) as dst:
    dst.write(result)

This produced the desired output. Note that my original image has RGB values in [0, 255] rather than in [0, 1]; hence, the dtype value of profile was uint8 in the single-band GTiffs and the calculation resulted in floats, which is why I added the line profile["dtype"] = "float64" to have the correct dtype to write the output (maybe useful for some reader).

Since this part appears to work, I assume you lost the geoinformation when you split the original image into R, G, B, alpha bands? This should be no problem, since you could simply obtain the profile from the original GTiff instead of the red band. They should only differ in the number of bands, which can be fixed easily:

with rasterio.open(PATH_TO_ORIGINAL_GTIFF) as f:
    profile = f.profile
# set the number of bands in the extracted information to 1
profile["count"] = 1

Complete:

import rasterio
with rasterio.open('path/to/Original_4_channel.tif') as f:
    profile = f.profile
# set the number of bands to 1
profile["count"] = 1

# [INSERT WHAT EVER CODE YOU USED TO GET THE SINGLE R, G AND B CHANNELS HERE]
# e.g. load the image as img and run
# (red, green, blue, alpha) = np.transpose(img, axes = (2,0,1))

with rasterio.open('path/to/Red_channel.tif') as f:
    red = f.read()

with rasterio.open('path/to/Blue_channel.tif') as f:
    blue = f.read()

result = ((red**2)+(blue**2))/(blue) 

with rasterio.open('path/to/Output.tif', 'w', **profile) as dst:
    dst.write(result)

Edit

It was not part of your original question. Nevertheless, I included a way to get the separate raster channels before doing the calculations. This is a complete script that just takes the directory and file name of the original 4-channel GTiff image.

wd = "/path/to/your/GTiff_file"# just the folder without the filename
import rasterio, os
import numpy as np
img_path = os.path.join(wd, "INPUT_FILE_NAME.tif")# filename goes in here
out_path = os.path.join(wd, "OUTPUT_FILE_NAME.tif")

with rasterio.open(img_path) as f:
    profile = f.profile
    R = f.read(1)
    G = f.read(2)
    B = f.read(3)
    a = f.read(4)

profile["count"] = 1
profile["dtype"] = "float64"

result = ((R**2)+(B**2))/(B)
result = np.expand_dims(result, axis=0)

with rasterio.open(out_path, 'w', **profile) as dst:
    dst.write(result)
Related Question