#!/usr/bin/env python

import gdal
import osr
import sys
import osr

# =============================================================================
def Usage():
    print 'Usage: create_mosaic.py <coordinate_sys> <mosaic> <raster> ...'
    sys.exit( 1 )

# =============================================================================
# Program mainline.
#

if len(sys.argv) < 4:
    Usage()

srs_arg = sys.argv[1]

mosaic = sys.argv[2]
files = sys.argv[3:]

# =============================================================================
# Instantiate the coordinate system.

target_srs = osr.SpatialReference()
if target_srs.SetFromUserInput( srs_arg ) != 0:
    print 'Failed to create SRS.'
    sys.exit( 1 )

print 'Generating Mosaic in Coordinate System:'
print target_srs.ExportToPrettyWkt()

# =============================================================================
# Process all files.

targ_xmin = None
targ_ymin = None
targ_xmax = None
targ_ymax = None

for filename in files:

    fd = gdal.Open( filename )

    # Read projection, and create transformer to target SRS. 
    file_srs_wkt = fd.GetProjection()
    file_srs = osr.SpatialReference()
    file_srs.ImportFromWkt( file_srs_wkt )
    
    ct = osr.CoordinateTransformation( file_srs, target_srs )
    
    # Setup bounds in source srs.
    gt = fd.GetGeoTransform()
    
    ulx = gt[0]
    uly = gt[3]
    
    lrx = gt[0] + fd.RasterXSize * gt[1] + fd.RasterYSize * gt[2]
    lry = gt[3] + fd.RasterXSize * gt[4] + fd.RasterYSize * gt[5]
    
    corners = [ (ulx, uly), (ulx, lry), (lrx, uly), (lrx, lry) ]
    
    for corner in corners:
        (x,y,z) = ct.TransformPoint( corner[0], corner[1] )
        if targ_xmin is None:
            targ_xmin = x
            targ_xmax = x
            targ_ymin = y
            targ_ymax = y
        else:
            targ_xmin = min(targ_xmin,x)
            targ_xmax = max(targ_xmax,x)
            targ_ymin = min(targ_ymin,y)
            targ_ymax = max(targ_ymax,y)

if targ_xmin is None:
    print 'Failed to get extents, stopping.'
    sys.exit( 1 )

# =============================================================================
# What is the target file size and geotransform?
#
# For now we will hardcode 60m x 60m pixels.

x_pixel_size = 60.0
y_pixel_size = 60.0

xsize = (targ_xmax - targ_xmin) / x_pixel_size + 1
ysize = (targ_ymax - targ_ymin) / y_pixel_size + 1

targ_gt = (targ_xmin, x_pixel_size, 0.0,
           targ_ymax, 0.0, - y_pixel_size)

# =============================================================================
# Create target RGB "mosaic" file.  

driver = gdal.GetDriverByName( 'GTiff' )

fd = driver.Create( mosaic, xsize, ysize, 3, gdal.GDT_Byte, [ 'TILED=YES' ] )

fd.SetProjection( target_srs.ExportToWkt() )
fd.SetGeoTransform( targ_gt )