Search code examples
geospatialgdalrasterio

geotiff exported all black with rasterio

I'm developing a workflow for satellite imagery processing with Python, using rasterio (1.3.6) and earthpy libraries.

The idea is to generate some raster algebra, make a composition (sabinsRatio), and then export this composed raster as one file, where each band has one ratio.

I can see the output in the Python env, but when I try to see it in QGIS, it only displays a black square. I have looked at several questions from Stackoverflow and the documentation, but still couldn't figure it out.


#file path
multi_bands = glob.glob('.././GIS/landsat_8/LC08_L2SP_090079_20230225_20230301_02_T1/LC08_L2SP_090079_20230225_20230301_02_T1_SR_*B[2:3:4:5:6:7].tif')

multi_bands.sort()

img_list = []

#reading bands
for img in multi_bands:
    with rio.open(img, 'r') as img_file:
        img_list.append(img_file.read(1))

#multiband array    
arr_st = np.stack(img_list)

# visualize all bands
ep.plot_bands(arr_st,
              cbar=False,
              cmap = 'gist_earth')
plt.show()

#avoid numpy error of zero division
np.seterr(divide='ignore', invalid='ignore')

#perform band ratios
ironOxideR = arr_st[0]/arr_st[2]
ClayHyR = arr_st[4]/arr_st[5]
ferrousR = arr_st[4]/arr_st[3]

# Create the RGB composition
sabinsRatio = np.stack((ironOxideR, ClayHyR, ferrousR))

#sabinsRatio plot
ep.plot_rgb(sabinsRatio,
            rgb=(0,1,2),
            figsize=(10, 10),
            stretch = True,
            title = 'RGB composition of Iron Oxide, Clay Hydration, and Ferrous')
plt.show()

#trying to export

#exporting sabinsRatio

# Update the metadata
out_meta = {"driver": "GTiff",
            "height": sabinsRatio.shape[1],
            "width": sabinsRatio.shape[2],
            "crs": "EPSG:32756",
            "count":3,
            "dtype":rio.float32
            }
           


# Write the mosaic raster to disk
with rio.open('./output/sabinRation2.tif', "w", **out_meta) as dest:
    dest.write(sabinsRatio, [1,2,3])

1. The plot of all bands enter image description here

2. Sabins Ratio plot enter image description here

3. QGIS Visualization: if you see, there's a small red dot in the upper left corner of the black square. That's my area of interest. enter image description here

How could I fix this?

Solution

  • I found the solution with the aid of a friend. In total, two changes were made:

    First

    
for img in multi_bands:
    with rio.open(os.path.join(path,img), 'r') as img_file:
        img_list.append(img_file.read(1))
        out_meta = img_file.meta.copy()

    Second

    
# Update the metadata
out_meta.update({"driver": "GTiff",
            "height": sabinsRatio.shape[1],
            "width": sabinsRatio.shape[2],
            "count":3,
            "dtype":rio.float32
            })

    Finally

    it was possible to open on QGIS.