4

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.

6

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))
  • 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 '15 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 '15 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 '15 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 '15 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 '15 at 21:03

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.