Overlaying Sentinel 2 13 Band Data with Classification TIF Data in Python

Importing Packages

We will utilize 4 key Python libraries:

import rasterio
from rasterio.plot import show
import earthpy.spatial as es
import matplotlib.pyplot as plt

• Rasterio: reading and manipulating GeoTIFF/raster data
• Earthpy: satellite and geospatial processing tools
• Matplotlib: visualization and plotting capabilities

Loading Sentinel 2 Imagery

We will demonstrate with a sample Sentinel 2 L1C tile covering part of Belgium:

s2_path = 's2_tile.tif'
src = rasterio.open(s2_path)
s2 = src.read()
print(f'Sentinel bands: {s2.shape}')

This loads the 13 spectral bands as a NumPy array with dimensions: (13, 10980, 10980)”

Load Classification TIF Data

Next, we open the classification raster for the same region:

ml_path = 'crop_map.tif' 
ml_src = rasterio.open(ml_path)
ml_arr = ml_src.read()
print(f'Classification array shape: {ml_arr.shape}')

Our crop_map is a single-channel TIF with categorical values.

Align GeoTIFF Extents

Now we clip the Sentinel 2 coverage to match our classified TIF layer exactly:

aligned_s2_arr, ml_transform = es.crop_image(s2, ml_src)

This ensures the same spatial extent and band counts when overlaying.

Map RGB Bands

Next, we define which Sentinel bands map to RGB channels for eventual visualization:

red = aligned_s2_arr[3].squeeze() 
green = aligned_s2_arr[2].squeeze()
blue = aligned_s2_arr[1].squeeze()

Band 4 covers red wavelengths, band 3 covers green, and band 2 captures the blue channel.

Define Color Palette

Then we create a Categorical colormap to colorize pixels based on classified TIFF categories:

cmap = plt.colormaps['Pastel2'] 
colors = [ cmap(i) for i in range(cmap.N) ]

Apply Colors

Here we map RGB channels and colorize where classification exists:

vis = ml_arr.copy()
vis[ml_arr==0] = (red[ml_arr==0], green[ml_arr==0],  
                  blue[ml_arr==0]) 

vis[ml_arr==1] = colors[1] 

vis[ml_arr==2] = colors[2]

Areas where no classification label is assigned retain true color from S2. Other regions get colorized based on their coding value to reveal segmentation.

Visualize Output

Finally, plot using Matplotlib:

fig, ax = plt.subplots(figsize=(12, 12))  
ax.imshow(vis) 
ax.set_title('Sentinel 2 overlaid with Crop Map')
ax.set_axis_off() 
plt.show()

We should now see an integrated raster with a true color background from Sentinel 2 images enriched by classified TIF overlays!

Full Script

Here is the entire overlay code workflow:

import rasterio
from rasterio.plot import show
import earthpy.spatial as es
import matplotlib.pyplot as plt

s2_path = 's2_tile.tif’  
src = rasterio.open(s2_path)
s2 = src.read() 

ml_path = 'crop_map.tif’   
ml_src = rasterio.open(ml_path)  
ml_arr = ml_src.read()   

aligned_s2_arr, ml_transform = es.crop_image(s2, ml_src)
   
red = aligned_s2_arr[3].squeeze()  
green = aligned_s2_arr[2].squeeze()
blue = aligned_s2_arr[1].squeeze()   

cmap = plt.colormaps['Pastel2']  
colors = [ cmap(i) for i in range(cmap.N) ]
   
vis = ml_arr.copy()
vis[ml_arr==0] = (red[ml_arr==0], green[ml_arr==0], blue[ml_arr==0])
   
vis[ml_arr==1] = colors[1]  
vis[ml_arr==2] = colors[2]   

fig, ax = plt.subplots(figsize=(12, 12))   
ax.imshow(vis)
ax.set_title('Sentinel 2 overlaid with Crop Map')
ax.set_axis_off()
plt.show()

Automating this Earth observation data fusion process enables the rapid generation of insightful landcover change visualizations!

Key Takeaways

To recap the main steps covered:

• Open Sentinel 2 tiles and classified TIFs with Rasterio
• Align geo extents between both datasets
• Define RGB channel mappings
• Create a color palette from classification codes
• Overlay truecolor background with categorical overlays
• Visualize integrated raster mosaic

This workflow can extend towards video generation by iterating across time series imagery with class raster animations.

We have only touched the surface of deep multivariate geospatial analytics possible from these types of multispectral satellite and vector data integrations within Python’s powerful computational geospatial libraries.

But this initial tutorial should provide a template to get you started overlaying additional datasets for all kinds of insightful Earth analytics!