Write raster file using rasterio in Python?

Question

I am trying to follow a tutorial for computing NDVI (Normalized Difference Vegetation Index) through the rasterio package in python, however, I am unsure how to finish the task by actually creating the raster .tif file itself. I am following this tutorial here: https://automating-gis-processes.github.io/CSC18/lessons/L6/raster-calculations.html

from Henrikki Tenkanen.

I follow the tutorial all the way until the end and then am left with ndvi which I can see is an array, and I can find its mean and also plot/map it like shown in the tutorial. However, how do I actually convert this object/array to a raster .tif file that I can actually view in QGIS? It seems like that step is missing from the tutorial.

Just looking through rasterio documentation, it seems like there is some need to use a dst.write() argument, but I am just too new at this to really understand how to actually use this process. It looks something like:

with rasterio.open('example.tif', 'w', **profile) as dst:
     dst.write(array.astype(rasterio.uint8), 1)

Here is the dst.write() documentation I got this from: https://rasterio.readthedocs.io/en/latest/topics/writing.html

I cannot figure out how to send ndvi to a .tif raster with the correct syntax. How can I complete this tutorial so that an actual ndvi raster .tif file is created? I am confused by what the kwargs are and how to use them.

Answer 1

I would recommend using @radouxju approach to calculating NDVI in this answer. Here is an untested example based on the link you provided and radouxju's NDVI approach:

import os
import rasterio
import numpy as np

# Filepath
dataset = r'C:\HY-DATA\HENTENKA\CSC\Data\Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif'

# Outfile path
outpath = r'C:\path\to\your\dir'

# Open multiband raster
raster = rasterio.open(dataset)

# Calc NDVI
red = raster.read(4)
nir = raster.read(5)

ndvi = np.zeros(red.shape, dtype=rasterio.float32)
ndvi = (nir.astype(float)-red.astype(float))/(nir+red)

# Write to TIFF
kwargs = red.meta
kwargs.update(
    dtype=rasterio.float32,
    count=1,
    compress='lzw')

with rasterio.open(os.path.join(outpath, 'ndvi.tif'), 'w', **kwargs) as dst:
    dst.write_band(1, ndvi.astype(rasterio.float32))

Answer 2

# Look At this code:

import rasterio
import numpy as np
import matplotlib.pyplot as plt

# open the bands:

B4 = rasterio.open(inputpath_name)
B5 = rasterio.open(inputpath_name)

# First calculate your NDVI with this function:

def ndvi(red, nir):
    red  = B4.read(1)
    nir  = B5.read(1)
    red = red.astype('float64')
    nir = nir.astype('float64')
    ndvi = (nir - red)/(nir + red)
    return ndvi


ndvi = ndvi(B4, B5) # Choose the bands inside the function and run

# plot your ndvi:

plt.figure(figsize=(9,9))
plt.imshow(ndvi, cmap='RdYlGn')
plt.colorbar()
plt.title('NDVI - Landsat 8 OLI')
plt.xlabel('Column #')
plt.ylabel('Row #')

# Edit your metadata like a dictionary in python:

out_meta = B4.meta.copy()

out_meta.update({'driver':'GTiff',
                 'width':B4.shape[1],
                 'height':B4.shape[0],
                 'count':1,
                 'dtype':'float64',
                 'crs':B4.crs, 
                 'transform':B4.transform,
                 'nodata':0})

# then:

# you have two options to save the index in format Geotiff using rasterio:


with rasterio.open(fp=r'ndvi.tif', #outputpath_name
             mode='w',**out_meta) as dst:
             dst.write(ndvi, 1) # the numer one is the number of bands

with rasterio.open(fp=r'ndvi.tif', # outputpath_name
             mode='w',**out_meta) as dst:
             dst.write_band(1,ndvi) # the numer one is the number of bands