1

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.

1

2 Answers 2

0

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))
0

# 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.