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.

• Kartograaf: have you seen this thread? -> gis.stackexchange.com/questions/138914/…

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