from math import ceil
import os, sys

from owslib.wms import WebMapService
from pyproj import Proj, transform
from PIL import Image

wms = WebMapService(sys.argv[1], version='1.1.1')

if len(sys.argv) < 3:
    first_layer = None
    for layer in wms.contents:
        if first_layer is None:
            first_layer = layer
        print layer
    print wms[first_layer].boundingBox
    print wms[first_layer].boundingBoxWGS84
    sys.exit()
else:
    layers = sys.argv[8:]
    bbox = [float(b) for b in sys.argv[2:6]]
    coord_sys = wms[layers[0]].boundingBox[4]
    wgs = Proj(init='epsg:4326')
    out_proj = Proj(init=coord_sys)
    result_bbox = transform(wgs, out_proj, bbox[0], bbox[1]) + transform(wgs, out_proj, bbox[2], bbox[3])

    image_size = [int(sys.argv[6]), int(sys.argv[7])]
    bbox_size = [abs(result_bbox[2] - result_bbox[0]), abs(result_bbox[3] - result_bbox[1])]
    bbox_ratio = bbox_size[0] / bbox_size[1]
    image_size_corrected = [
        int(max(image_size[0], image_size[1]*bbox_ratio)),
        int(max(image_size[1], image_size[0]/bbox_ratio))]
    print 'Correcting image size to %dx%d due to coordinates box ratio' % (
        image_size_corrected[0], image_size_corrected[1])

    image_size = [int(ceil(i/400.0))*400 for i in image_size_corrected]
    print 'Correcting image size to %dx%d due to tile size (400x400)' % (
        image_size[0], image_size[1])

    min_long = min(result_bbox[0], result_bbox[2])
    max_long = max(result_bbox[0], result_bbox[2])
    long_span = max_long - min_long
    min_lat = min(result_bbox[1], result_bbox[3])
    max_lat = max(result_bbox[1], result_bbox[3])
    lat_span = max_lat - min_lat
    tile_long_span = long_span * 400.0 / image_size_corrected[0]
    tile_lat_span = lat_span * 400.0 / image_size_corrected[1]

    tiles = []
    current_lat = min_lat
    current_long = min_long
    while current_lat < max_lat:
        current_row = []
        while current_long < max_long:
            current_row.append([
                current_long, current_lat,
                current_long+tile_long_span, current_lat+tile_lat_span])
            current_long += tile_long_span
        if len(current_row):
            tiles.append(current_row)
        current_lat += tile_lat_span
        current_long = min_long

    files = []
    x = 0
    y = 0
    for row in tiles:
        filerow = []
        for tile in row:
            print 'Fetching tile %dx%d: %s' % (x, y, tile)
            img = wms.getmap(layers=layers,
                             srs=coord_sys,
                             bbox=tile,
                             size=(400,400),
                             format='image/png',
                             transparent=True)
            filename = 'out-%d-%d.png' % (x, y)
            out = open(filename, 'wb')
            out.write(img.read())
            out.close()
            filerow.append(filename)
            y += 1
        files.append(filerow)
        x += 1
        y = 0

    strips = []
    r = 0
    for row in files:
        tiles = [Image.open(f) for f in row]
        row_image = Image.new('RGBA', (400*len(tiles), 400))
        offset = 0
        for tile in tiles:
            row_image.paste(tile, (offset, 0))
            offset += 400
        strip_filename = 'out-%d.png' % (r)
        row_image.save(strip_filename)
        strips.append(strip_filename)
        for f in row:
            os.remove(f)
        r += 1

    strips.reverse()
    strip_images = [Image.open(f) for f in strips]
    final_image = Image.new('RGBA', (strip_images[0].width, 400*len(strip_images)))
    offset = 0
    for strip_image in strip_images:
        final_image.paste(strip_image, (0, offset))
        offset += 400
    final_image.save('out.png')
    for s in strips:
        os.remove(s)