0

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
14
  • 1
    I have no experience using rasterio whatsoever, but my guess would be you only extract the values with your code and you would have to manually copy the origin and projection from the original image Oct 27, 2021 at 13:14
  • 1
    f.read() reads the raster data to a numpy array. Hence, geospatial information is lost at this point. Oct 27, 2021 at 13:26
  • 1
    Maybe this helps: medium.com/@mommermiscience/… Oct 27, 2021 at 13:27
  • 1
    If you have GDAL installed, you can use the code from this answer to get the geospatial information from the original file to your result gis.stackexchange.com/a/327239/150644 Oct 27, 2021 at 13:37
  • 1
    I just tried to reproduce the issue but I am not sure where np.transpose comes in. I split a GTiff into R, G, and B and used your code to execute the calculation and save the output. QGIS displayed the resulting GTiff right above the original one, so I assume the georeferences are correct. Did you use np.transpose to generate the separate bands from a 4-band tiff? Oct 28, 2021 at 11:06

1 Answer 1

1

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)
7
  • Please see my update on the initial post Oct 28, 2021 at 18:14
  • 1
    ValueError: axes don't match array means the image currently stored under the variable name f has a number of dimensions that doesn't match what the function expects, given the parameters you set. Maybe you could provide the image somehow or show how you load it and what properties it has? Oct 28, 2021 at 19:23
  • You are amazing!!!! The code works! One last question... in case I want to read the full path in the img_path, so the path + tif file, such as: "/path/to/your/GTiff_file/INPUT_FILE_NAME.tif" how I pass it to os.path.join(..). I tried this command: img_path = os.open("/path/to/your/GTiff_file/INPUT_FILE_NAME.tif", os.O_RDONLY) but it does not work.... Any ideas? Oct 29, 2021 at 11:00
  • 1
    If I understand correctly, you want to create img_path without using the wd variable? You could just replace wd with the directory path: img_path = os.path.join("/path/to/your/GTiff_file", "INPUT_FILE_NAME.tif") or you could set the variable directly as img_path = "/path/to/your/GTiff_file/INPUT_FILE_NAME.tif" Oct 29, 2021 at 11:10
  • 1
    The description of np.transpose says "Reverse or permute the axes of an array", so yes, you could probably change the order of bands. However, in the solution I gave there is no need for this function. Unfortunately, I just started learning Python 4 months ago and don't know every detail; I guess you will find better answers to such questions if you start a new question instead of commenting on this answer where few people will read it... I guess the people who know all the answers won't read this conversation Nov 7, 2021 at 19:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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