#!/usr/bin/env python

import glob
import string
import os
import time
import stat
import Image
import sys
import gdal
import gdalconst

chunked_dir = "data/square_chunked"
pixels_per_degree = 1350

# Read all of our chunks.
chunk_names = glob.glob('%s/*.tif' % (chunked_dir))

h_min = None
chunk_positions = []
for chunk_filename in chunk_names:
	chunk_name = os.path.basename(chunk_filename)[:-4]
	print chunk_name
	(junk, h_s, v_s) = chunk_name.split('_')
	h = int(h_s)
	v = int(v_s)

	chunk_positions.append((h, v))

	if h_min is None:
		h_min = h_max = h
		v_min = v_max = v
	else:
		h_min = min(h_min, h)
		h_max = max(h_max, h)
		v_min = min(v_min, v)
		v_max = max(v_max, v)

print h_min, h_max, v_min, v_max

# Open the first one to get a dimension.
print chunk_names[0]
#im = Image.open(chunk_names[0], 'r')
im = gdal.Open(chunk_names[0])
square_width = im.RasterXSize
square_height = im.RasterYSize
im = None

span_X = h_max - h_min + 1
span_Y = v_max - v_min + 1
print span_X, span_Y, span_X * span_Y

image_width = span_X * square_width
image_height = span_Y * square_height

print image_width, image_height

target_image = Image.new('RGB', (image_width, image_height))

def add_chunk(longitude, latitude):
	chunk_filename = chunked_dir + '/chunk_%d_%d.tif' % (longitude, latitude)
	print 'add %s' % (chunk_filename)

	x_offset = (longitude - h_min) * pixels_per_degree
	y_offset = (latitude - v_min) * pixels_per_degree
	print x_offset, y_offset

	chunk_image = gdal.Open(chunk_filename, gdalconst.GA_ReadOnly)

	for band in range(3):
		target_image.GetRasterBand(band + 1).WriteArray(
			chunk_image.GetRasterBand(band + 1).ReadAsArray(0, 0, square_width, square_height),
			x_offset, y_offset
		)

for (lon, lat) in chunk_positions:
	print lon, lat
	add_chunk(lon, lat)