Test script to flatten some spikes in some aster data.

flatten_spikes.py

from osgeo import gdal
import numpy as np
import os

ds = gdal.Open("UNIT_S60W070/ASTGTM2_S56W068_dem.tif")

band = ds.GetRasterBand(1)

data = band.ReadAsArray()

threshold = 1000.0

total = data.shape[0] * data.shape[1]
num_scaled = 0

for r in range(data.shape[0]):
	for c in range(data.shape[1]):
		h = data[r,c]
		if h > threshold:
			#h = h * scale
			# Actually flatten, don't scale the pixel.  This looks better.
			h = 0
			data[r,c] = h
			num_scaled += 1
			
print "Scaled %s of %s pixels" % (num_scaled, total)
# Create an output geotiff to hold the corrected data.
destination = "no_spikes.tif"
if os.path.exists(destination):
	os.unlink(destination)
out_ds = gdal.GetDriverByName("GTiff").Create(destination, band.XSize, band.YSize, 1, gdal.GDT_Int16, ['COMPRESS=LZW'])
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_ds.SetProjection(ds.GetProjection())             
out_ds.GetRasterBand(1).WriteArray(data, 0, 0)