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