#!/usr/bin/env python
###############################################################################
# $Id$
#
# Purpose:  Builds global VRT for GMTED2010 data.
# Author:   Even Rouault, <even dot rouault at mines dash paris dot org>
#
###############################################################################
# Copyright (c) 2011, Even Rouault, <even dot rouault at mines dash paris dot org>
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
###############################################################################

from osgeo import gdal
import os

def build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon):
    return '/vsicurl/http://igskmncngs506.cr.usgs.gov/gmted/Global_tiles_GMTED/%sdarcsec/%s/%s%03d/%02d%s%03d%s_20101117_gmted_%s%s.tif' % (res, typ, lon_hemisphere, abs_lon, abs_lat, lat_hemisphere, abs_lon, lon_hemisphere, typ, res)

def build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon):
    return 'gmted/%02d%s%03d%s_20101117_gmted_%s%s.vrt' % (abs_lat, lat_hemisphere, abs_lon, lon_hemisphere, typ, res)

def add_overviews(filename, ovr_name_list):
    f = open(filename, 'rt')
    lines = f.readlines()
    f.close()

    for line in lines:
        if line.find('<Overview') != -1:
            return

    f = open(filename, 'wt')
    for line in lines:
        if line.find('</VRTRasterBand>') != -1:
            for ovr_name in ovr_name_list:
                if ovr_name.find('/vsicurl/') == -1:
                    relative = 1
                else:
                    relative = 0
                f.write('    <Overview>\n')
                f.write('      <SourceFilename relativeToVRT="%d">%s</SourceFilename>\n' % (relative, ovr_name))
                f.write('      <SourceBand>1</SourceBand>\n')
                f.write('    </Overview>\n')
        f.write('%s' % line)
    f.close()

def build_vrt(vrt_name, tile_list):
    f = open('vrt_list.txt', 'wt')
    for filename in tile_list:
        f.write('%s\n' % filename)
    f.close()
    os.system('gdalbuildvrt -resolution highest -input_file_list vrt_list.txt %s' % vrt_name)
    os.unlink('vrt_list.txt')


try:
    os.stat('gmted')
except:
    os.mkdir('gmted')

gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS',  ".tif")

vrt_list075 = []
vrt_list150 = []
vrt_list300 = []

lon_step = 30
lat_step = 20

lat = 90 - lat_step
while lat >= -90:
    abs_lat = abs(lat)
    if lat < 0:
        lat_hemisphere = 'S'
    else:
        lat_hemisphere = 'N'

    lon = -180
    while lon < 180:
        abs_lon = abs(lon)
        if lon < 0:
            lon_hemisphere = 'W'
        else:
            lon_hemisphere = 'E'

        # Special case for southern most tiles : there's only data at resolution 300 and for mea type
        if lat == -90:
            typ = 'mea'
            res = '300'
        else:
            # We can use the following values for the type :
            # bln: Breakline Emphasis
            # dsc: Systematic Subsample
            # mea: Mean
            # med: Median
            # min: Minimum
            # max: Maximum
            # std: Standard Deviation
            typ = 'bln'
            res = '075'

        src_name = build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
        target_name = build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)

        try:
            os.stat(target_name)
            vrt_list075.append(target_name)
        except:
            print('Generating %s' % target_name)
            src_ds = gdal.Open(src_name)
            if src_ds is not None:
                vrt_ds = gdal.GetDriverByName('VRT').CreateCopy(target_name, src_ds)
                vrt_ds = None
                vrt_list075.append(target_name)

        if target_name in vrt_list075 and res == '075':
            ovr_name_list = []
            for ovr_res in ['150', '300']:
                ovr_name = build_src_filename(ovr_res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
                ovr_name_list.append(ovr_name)

            add_overviews(target_name, ovr_name_list)

            res = '150'
            src_name = build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
            target_name = build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)

            try:
                os.stat(target_name)
                vrt_list150.append(target_name)
            except:
                print('Generating %s' % target_name)
                src_ds = gdal.Open(src_name)
                if src_ds is not None:
                    vrt_ds = gdal.GetDriverByName('VRT').CreateCopy(target_name, src_ds)
                    vrt_ds = None
                    vrt_list150.append(target_name)

            res = '300'
            src_name = build_src_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)
            target_name = build_target_filename(res, typ, lat_hemisphere, abs_lat, lon_hemisphere, abs_lon)

            try:
                os.stat(target_name)
                vrt_list300.append(target_name)
            except:
                print('Generating %s' % target_name)
                src_ds = gdal.Open(src_name)
                if src_ds is not None:
                    vrt_ds = gdal.GetDriverByName('VRT').CreateCopy(target_name, src_ds)
                    vrt_ds = None
                    vrt_list300.append(target_name)

        elif target_name in vrt_list075 and res == '300':
            vrt_list150.append(target_name)
            vrt_list300.append(target_name)

        lon += lon_step

    lat -= lat_step

build_vrt('gmted/all075.vrt', vrt_list075)
build_vrt('gmted/all150.vrt', vrt_list150)
build_vrt('gmted/all300.vrt', vrt_list300)

add_overviews('gmted/all075.vrt', ['all150.vrt', 'all300.vrt'])
add_overviews('gmted/all150.vrt', ['all300.vrt'])