#!/usr/bin/env python

import sys
import gdal
import gdalconst

chunk_size_X = chunk_size_Y = 100000

input_filename = sys.argv[1]

print "input_filename=%s" % (input_filename)

input_image = gdal.Open(input_filename, gdalconst.GA_ReadOnly)

# (-777482.09852564591, 64.0, 0.0, -164027.04259051729, 0.0, -64.0)
input_geotransform = input_image.GetGeoTransform()
print 'input_geotransform=%s ' % (str(input_geotransform))

#(ul_corner_X, pixel_size_X, skew, ul_corner_Y, pixel_size_Y, zero, pixel_size_Y) 

input_left_X = input_geotransform[0]
input_top_Y = input_geotransform[3]
pixel_size_X = input_geotransform[1]
pixel_size_Y = input_geotransform[5]

input_right_X = input_geotransform[0] + input_geotransform[1] * input_image.RasterXSize
input_lower_Y = input_geotransform[3] + input_geotransform[5] * input_image.RasterYSize

print 'LR corner=%f, %f' % (input_right_X, input_lower_Y)

# Start in the upper left corner (pixel 0,0).
AOI_X = input_left_X
AOI_Y = input_top_Y

while 1:
	# print 'location=%0.1f, %0.1f' % (AOI_X, AOI_Y),

	chunk_X = int(AOI_X / chunk_size_X) * chunk_size_X
	chunk_Y = int(AOI_Y / chunk_size_Y) * chunk_size_Y

	# print 'chunk location=%0.1f, %0.1f' % (chunk_X, chunk_Y),

	chunk_extent_X = chunk_X + chunk_size_X
	if chunk_extent_X > input_right_X:
		chunk_extent_X = input_right_X
		next_X = input_left_X
		next_Y = chunk_extent_Y	
	else:
		next_X = chunk_extent_X
		next_Y = AOI_Y

	chunk_extent_Y = chunk_Y - chunk_size_Y

	if chunk_extent_Y < input_lower_Y:
		chunk_extent_Y = input_lower_Y

	# print 'chunk_extent=%0.1f, %0.1f' % (chunk_extent_X, chunk_extent_Y),

	print AOI_X, AOI_Y, chunk_extent_X, chunk_extent_Y

	if chunk_extent_X >= input_right_X and chunk_extent_Y <= input_lower_Y:
		print
		print 'breaking'
		print chunk_extent_X, input_right_X, chunk_extent_Y, input_lower_Y
		break

	# Move to the next chunk.
	AOI_X = next_X
	AOI_Y = next_Y