7

I'm using Rasterio to read GeoTIFF files from Landsat 8 and calculate NDVI into a new GeoTIFF file.

My code looks like this:

import numpy
import rasterio
import subprocess

with rasterio.drivers(CPL_DEBUG=True):

    # Read raster bands directly to Numpy arrays.
    #
    dsRed = rasterio.open('downloads/LC81970212013122LGN01/LC81970212013122LGN01_B3.TIF')
    bandRed = dsRed.read_band(1)

    dsNir = rasterio.open('downloads/LC81970212013122LGN01/LC81970212013122LGN01_B5.TIF')
    bandNir = dsNir.read_band(1)

    ndvi = numpy.zeros(dsRed.shape, dtype=rasterio.uint16)

    ndvi_upper = bandNir + bandRed
    ndvi_lower = bandNir - bandRed

    ndvi = ndvi_lower / ndvi_upper

    kwargs = dsRed.meta
    kwargs.update(
        dtype=rasterio.uint16,
        count=1,
        compress='lzw')

    with rasterio.open('example-total.tif', 'w', **kwargs) as dst:
        dst.write_band(1, ndvi.astype(rasterio.uint16))


# Dump out gdalinfo's report card and open the image.
info = subprocess.check_output(['gdalinfo', '-stats', 'example-total.tif'])
print(info)

subprocess.call(['open', 'example-total.tif'])

But it isn't really producing the result I was hoping for.

I'm probably getting the data types wrong, I'm not too familiar with Python and NumPy.

The resulting image is all black. If I multiple the result array by e.g. 1000 it becomes visible, but it isn't the correct values.

1 Answer 1

10

First I would use bands 4(red) and 5(nir) for Landsat 8 according to the description of the OLI instrument, and 3(red) and 4(NIR) for the Landsat TM and ETM.

Second, you define an output in dtype=rasterio.uint16, but NDVI should be a float (between -1 and 1). You should either initialize your raster as dtype=rasterio.float32 , or multiply your values by 1000 in your code, ans store as int16.

ndvi = numpy.zeros(dsRed.shape, dtype=rasterio.float32)

ndvi = (bandNIR.astype(float)-bandRed.astype(float))/(bandNIR+bandRed)

kwargs = dsRed.meta
kwargs.update(
    dtype=rasterio.float32,
    count=1,
    compress='lzw')

with rasterio.open('example-total.tif', 'w', **kwargs) as dst:
    dst.write_band(1, ndvi.astype(rasterio.float32))
12
  • Works better using band 4 (I thought band 3 was red). Instead of being all black, all countryside is now dark. The same dark color all over. Is it supposed to be lower / upper, or is it upper / lower? Changing it the result seems more correct.
    – MartinHN
    Mar 15, 2015 at 19:56
  • it has to be with the dif on top. Note that dark does not necessarily mean that you have a problem : it could only be a DISPLAY issue.
    – radouxju
    Mar 15, 2015 at 20:18
  • @MartinHN Landsat 8 band 4 is red, Landsat 4-7 (TM and ETM+) band 3 is red - landsat.usgs.gov/band_designations_landsat_satellites.php
    – user2856
    Mar 15, 2015 at 20:28
  • @radouxju It probably is a display issue. For display it would also make more sense to visualise similar to this: upload.wikimedia.org/wikipedia/commons/9/94/NDVI_062003.png
    – MartinHN
    Mar 15, 2015 at 20:53
  • Although, looking at the stats of the resulting file it seems weird that none of the values are below zero? Metadata: STATISTICS_MAXIMUM=65000 STATISTICS_MEAN=1321.0759915166 STATISTICS_MINIMUM=0 STATISTICS_STDDEV=2195.2510197044
    – MartinHN
    Mar 15, 2015 at 21:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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