In [232]:
from __future__ import print_function, division
from datetime import datetime as dt
import fnmatch
import os
from osgeo import gdal
from osgeo import gdal_array
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Locate example dataset from: https://github.com/ceholden/landsat_stack.git
We need the names and the ordinal dates of each stacked image.
In [2]:
location = '/home/ceholden/Documents/landsat_stack/p022r049/images'
# Pattern for stack image directories
stack_pattern = 'L*'
# Pattern for stack images
image_pattern = 'L*stack'
In [30]:
dates = []
images = []
for dirname in fnmatch.filter(os.listdir(location), stack_pattern):
# Find date as datetime object
stack_date = dt.strptime(dirname[9:16], '%Y%j').toordinal()
# Prepend the location to dirname
dirname = os.path.join(location, dirname)
# Find the stack image
stack = fnmatch.filter(os.listdir(dirname), image_pattern)
if len(stack) != 1:
raise Exception('More than one or no stack images found in {d}'.format(d=dirname))
images.append(os.path.join(dirname, stack[0]))
dates.append(stack_date)
# Store as NumPy arrays
images = np.array(images)
dates = np.array(dates)
# Sort images by date
sort_index = np.argsort(dates)
dates = dates[sort_index]
images = images[sort_index]
print(dates.size, images.size)
In [31]:
date_min = dt.strptime('2000-01-01', '%Y-%m-%d').toordinal()
date_max = dt.strptime('2000-12-31', '%Y-%m-%d').toordinal()
Subset:
In [32]:
date_index = (dates <= date_max) & (dates >= date_min)
dates = dates[date_index]
images = images[date_index]
print(dates.size, images.size)
In [248]:
# First, identify the image attributes
example_ds = gdal.Open(images[0], gdal.GA_ReadOnly)
# Size
nrow = example_ds.RasterYSize
ncol = example_ds.RasterXSize
nband = example_ds.RasterCount
# Data type
dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
example_ds.GetRasterBand(1).DataType)
gdal_dtype = example_ds.GetRasterBand(1).DataType
# Save projection & geo-transform for output
proj = example_ds.GetProjection()
geo_transform = example_ds.GetGeoTransform()
example_ds = None
In [253]:
# Now, init storage array and read in
stack = np.zeros((images.size, nrow, ncol, nband))
for i, image in enumerate(images):
ds = gdal.Open(image, gdal.GA_ReadOnly)
for b in range(nband):
band = ds.GetRasterBand(b + 1)
stack[i, :, :, b] = band.ReadAsArray()
ds = None
In [238]:
def get_ndvi(stack):
return ((stack[0, :, :, 3] - stack[0, :, :, 2]) /
(stack[0, :, :, 3] + stack[0, :, :, 2]).astype(np.float))
In [241]:
ndvi = get_ndvi(stack)
plt.subplot(121)
plt.imshow(ndvi, cmap=plt.cm.gray)
plt.subplot(122)
hist = plt.hist(ndvi)
Zhe Zhu's compositing algorithm picks pixels from the various dates of images that have the maximum ratio of Band 4 (NIR) to Band 1 (blue)
In [234]:
max_vi = np.argmax(stack[:, :, :, 3] / stack[:, :, :, 1], axis=0)
In [235]:
k, j = np.meshgrid(np.arange(nrow), np.arange(ncol))
In [244]:
comp = stack[max_vi, j, k]
print(comp.shape)
In [243]:
ndvi = get_ndvi(comp[np.newaxis, :, :, :])
plt.subplot(121)
plt.imshow(ndvi, cmap=plt.cm.gray)
plt.subplot(122)
hist = plt.hist(ndvi)
In [254]:
# Form filename
fname = 'comp_' + str(dt.fromordinal(date_min).date()) + '2' + str(dt.fromordinal(date_max).date()) + '.gtif'
output = os.path.join(location, fname)
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output,
ncol, nrow, nband,
gdal_dtype)
for b in range(nband):
out_ds.GetRasterBand(b + 1).WriteArray(comp[:, :, b])
out_ds.SetProjection(proj)
out_ds.SetGeoTransform(geo_transform)
out_ds = None