Visualising geospatial .tiff images with Rasterio

Question

I am trying to visualise a .tiff image in Jupiter Notebook using Rasterio. I am a Junior Data Scientist for an AgriTech company and we just got access to 8 data layers (NDVI etc.) for two farms in .tiff format.

Here is the metadata for one image:

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -125125.0, 'width': 72, 'height': 87, 'count': 1, 'crs': CRS.from_epsg(32734), 'transform': Affine(20.0, 0.0, 364480.0,
       0.0, -20.0, 6292100.0), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'lzw', 'interleave': 'band'}

When I run the following:

ax = plt.figure(figsize=(15,10))

pic = rasterio.open('/content20180109_biow_Meerlust.tif','r',driver='GTiff',width=72,
      height=87,count=1, nodata=-125125.0)

show(pic,with_bounds=False)

I get a very pixelated image:

How do I visualise the image without it being so pixelated? My knowledge of the array adjustments behind these .tiff images is limited as I just started in the Agronomics field. Open to any suggestions.

My aim is to create a web app using Streamlit where I can overlay these images and create a short video of how the layers change over time.

Answer 1

here are a couple solutions that might help visualize multiple-band rasters with clarity. In both examples, raster is a rasterio.DatasetReader with multiple bands (indexed at 1).

1. Single Image

To view all layers in a single 2D plane, the bands have to be concatenated:

import numpy as np
import matplotlib.pyplot as plt

bands = []
for i in range(raster.count):
    bands.append(raster.read(i+1, out_dtype=np.intc))

plt.title("Full Color Raster")
plt.imshow(np.array(bands))
plt.show()

Unfortunately, because of limitations from pyplot's imshow() function, this method only works with a few layers (traditionally RGB). Feel free to use datatypes other than np.intc.

2. Visualize Layers Separately

The earthpy.plot module has several clean options for visualizing raster layers, including the convenient plot_bands():

import numpy as np
import earthpy.plot as ep

bands = []

# Read the raster's bands to an array
for i in range(raster.count):
    bands.append(raster.read((i+1), out_dtype=raster.dtypes[i]))

# Convert to an iterable np.ndarray and plot in a 3-column grid
ep.plot_bands(np.array(bands), cols=3) 

Really hope this helps! This is my first Stack Overflow response, so let me know if there's anything critical I omitted.