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: https://gis.stackexchange.com/questions/407821/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