create_dataset.py

import math
import glob
import os
import sys
from osgeo import gdal
from osgeo import gdalconst
from qgis.utils import iface
from qgis.core import QgsProject
from qgis.core import QgsApplication
from pathlib import Path
import shutil
from time import sleep
import multiprocessing
from multiprocessing import Pool
import concurrent.futures
import subprocess
import numpy
from itertools import repeat
from functools import partial
from datetime import datetime
from operator import itemgetter
import threading
from ftplib import FTP_TLS
from ftplib import FTP
from lxml import etree
from difflib import SequenceMatcher
import json
import io
import traceback

try:
    from PIL import Image
    import numpy
    import osgeo.gdal_array as gdalarray
    numpy_available = True
except ImportError:
    # 'antialias' resampling is not available
    numpy_available = False


CATEGORY = 'CreateDataset'

# enter one or more directory where your source files are
# ex. ['C:/Users/david/Documents/Minecraft/Source','C:/Users/david/Documents/Minecraft/Source2']
sources = ['C:/Users/david/Documents/Minecraft/Source'] 
# enter directory where you want to save the generated dataset/datasets. 
#If you wish to upload the dataset via ftp, change it to the path on your server (ex. Dataset/Tiled)
#Set it to None if you want to upload the dataset to root folder of your FTP/SFTP user
output = 'C:/Users/david/Documents/Minecraft/Tiled'
zoom = 15 #enter your zoom level
resampling_algorithm = 'near' # use near for most accurate color. Look at the resampling algorithm's comparison image on the wiki for other algorithms
manual_nodata_value = None #Leave at None to use the defined NODATA value of the source file, or set it to value, if your source files don't have NODATA defined
convert_feet_to_meters = False #Set to True, if your dataset heights are in feet

ftp_upload = False # Set to True for upload to FTP server
ftp_one_file = False #Set to True to upload one zip file (RenderedDataset.zip) to FTP server
ftp_s = False # Set to True, if you want to upload to a FTPS server
ftp_upload_url = '' # FTP url to upload to (Only IP address or domain, ex 192.168.0.26)
ftp_upload_port = 21 # FTP port, ex. 2121. Must be set
ftp_user = None # Leave at None for anonymous login, else set to user name, ex. 'davix'
ftp_password = None # Leave at None for anonymous login, else set to user password, ex. 'password'

cleanup = True # Set to False if you don't wish to delete VRT file and supporting files once the script completes. It will still do a cleanup, if you run the script again
thread_count = None #Set to the number of threads you want to use. Preferably don't use all threads at once. Leave at at None to use all threads

localThread = threading.local()

single_source_mode = None

try:
    import zipfile
except ImportError:
    ftpOnefile = False


# Modified version of scale_query_to_tile from gdal2tiles
#https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/scripts/gdal2tiles.py
class Tile:
  def __init__(self, ds, x, y, ulx, uly, lrx, lry, dsulx, dsxres, dsuly, dsyres, querysize=0):
        #dsulx, dsxres, dsxskew, dsuly, dsyskew, dsyres = ds.GetGeoTransform()
        rx = int((ulx - dsulx) / dsxres + 0.001)
        ry = int((uly - dsuly) / dsyres + 0.001)
        rxsize = max(1, int((lrx - ulx) / dsxres + 0.5))
        rysize = max(1, int((lry - uly) / dsyres + 0.5))

        if not querysize:
            wxsize, wysize = rxsize, rysize
        else:
            wxsize, wysize = querysize, querysize

        wx = 0
        if rx < 0:
            rxshift = abs(rx)
            wx = int(wxsize * (float(rxshift) / rxsize))
            wxsize = wxsize - wx
            rxsize = rxsize - int(rxsize * (float(rxshift) / rxsize))
            rx = 0
        if rx + rxsize > ds.RasterXSize:
            wxsize = int(wxsize * (float(ds.RasterXSize - rx) / rxsize))
            rxsize = ds.RasterXSize - rx

        wy = 0
        if ry < 0:
            ryshift = abs(ry)
            wy = int(wysize * (float(ryshift) / rysize))
            wysize = wysize - wy
            rysize = rysize - int(rysize * (float(ryshift) / rysize))
            ry = 0
        if ry + rysize > ds.RasterYSize:
            wysize = int(wysize * (float(ds.RasterYSize - ry) / rysize))
            rysize = ds.RasterYSize - ry

        self.x = x
        self.y = y
        self.rx = rx
        self.ry = ry
        self.rxsize = rxsize
        self.rysize = rysize
        self.wx = wx
        self.wy = wy
        self.wxsize = wxsize
        self.wysize = wysize
        self.querysize = querysize

class Job:
    def __init__(self, zoom_folder, zoom_level, input_file, bandsCount, resampling, output_path, ftpUpload, ftpOnefile, ftpS, ftpUrl, ftpPort, ftpUser, ftpPassword, datasetName):
        self.zoom_folder = zoom_folder
        self.zoom_level = zoom_level
        self.input_file = input_file
        self.bandsCount = bandsCount
        self.resampling = resampling
        self.outputPath = output_path
        self.ftpUpload = ftpUpload
        self.ftpOnefile = ftpOnefile
        self.ftpS = ftpS
        self.ftpUrl = ftpUrl
        self.ftpPort = ftpPort
        self.ftpUser = ftpUser
        self.ftpPassword = ftpPassword
        self.datasetName = datasetName
    
class TaskData:
    def __init__(self, sourceIndex, sourcesArray, timeStartedFirst):
        self.sourceIndex = sourceIndex
        self.sourcesArray = sourcesArray
        self.timeStartedFirst = timeStartedFirst

    def getDatasetName(self):
        return self.sourcesArray[self.sourceIndex].datasetName
    
    def getSourcePath(self):
        return self.sourcesArray[self.sourceIndex].sourcePath
    
    def setTilesCount(self, tiles_count):
        self.sourcesArray[self.sourceIndex].sourceTilesCount = tiles_count

    def getTilesCount(self):
        return self.sourcesArray[self.sourceIndex].sourceTilesCount
    
    def setTimeStarted(self, time_started):
        self.sourcesArray[self.sourceIndex].timeStarted = time_started

    def getTimeStarted(self):
        return self.sourcesArray[self.sourceIndex].timeStarted
    
    def setSourceBounds(self, wgs_minx, wgs_minz, wgs_maxx, wgs_maxz):
        self.sourcesArray[self.sourceIndex].sourceBounds = [wgs_minx, wgs_minz, wgs_maxx, wgs_maxz]

    def setSourceIndex(self, source_index):
        if source_index < len(self.sourcesArray):
            self.sourceIndex = source_index
            return True
        else:
            return False

    def genTaskName(self):
        return 'Create {dataset_name} elevation dataset ({index}/{length})'.format(dataset_name=self.getDatasetName(), index=self.sourceIndex + 1, length=len(self.sourcesArray))

class SourceDataset:
    def __init__(self, source_path):
        self.sourcePath = source_path
        self.datasetName = os.path.basename(source_path)
        self.sourceTilesCount = 0
        self.timeStarted = None
        self.sourceBounds = None

class ColorRamp:
    def __init__(self, altitude, cftm):
        v = None
        if cftm:
            v = altitude * 0.3048
        else:
            v = altitude
        v += 32768
        r = math.floor(v/256)
        g = math.floor(v % 256)
        b = math.floor((v - math.floor(v)) * 256)
        self.altitude = altitude
        self.red = r
        self.green = g
        self.blue = b

class Statistic:
    def __init__(self, minV, maxV):
        self.minV = minV
        self.maxV = maxV

class FTPArchiveProgress:
    def __init__(self, totalSize, task):
        self.totalSize = totalSize
        self.task = task
        self.sizeWritten = 0
        self.lastPr = 0

    def updateProgress(self, block):
        self.sizeWritten += 1024
        pr = max(0, min(int(((self.sizeWritten * 40) / self.totalSize) + 60), 100))
        if (self.lastPr != pr):
            self.lastPr = pr
            self.task.setProgress(pr)

def genTiles(task, taskData):
    taskData.setTimeStarted(datetime.now())
    try:
        cpu_count = 1
        if thread_count is None:
            cpu_count = multiprocessing.cpu_count()
        else:
            cpu_count = thread_count

        sources_length = len(taskData.sourcesArray)
        source_folder = taskData.getSourcePath()

        files = []
        for src in glob.iglob(source_folder + '/**', recursive=True):
            if (src.endswith(".tif") or src.endswith(".tiff")
            or src.endswith(".img") or src.endswith(".IMG")):
                src = src.replace("\\","/")
                fileinfo = QFileInfo(src)
                filename = fileinfo.completeBaseName()
                files.append([src, filename])

        if len(files) == 0:
            taskData.setTilesCount(0)
            return taskData

        if single_source_mode:
            QgsMessageLog.logMessage(
                        'Started creating dataset out of {count} files'.format(count=len(files)),
                        CATEGORY, Qgis.Info)
        else:
            QgsMessageLog.logMessage(
                        'Started processing {count} files in {dataset_name}({index}/{length}) to the combined output dataset'
                        .format(count=len(files), dataset_name=taskData.getDatasetName(), index=taskData.sourceIndex + 1, length=sources_length),
                        CATEGORY, Qgis.Info)
            

        org_vrt = os.path.join(source_folder, "OrgSource.vrt").replace("\\","/")
        if os.path.isfile(org_vrt):
            os.remove(org_vrt)
        
        org_files = []
        for file in files:
            try:
                tile_ds = gdal.Open(file[0], gdal.GA_ReadOnly)
                tile_info = gdal.Info(tile_ds, format='json')
                #Check if tiff has an positive NS resolution
                if tile_info["geoTransform"][5] > 0:
                    QgsMessageLog.logMessage('Reprojecting {tile_name} duo to positive NS resolution (vertically flipped image)'.format(tile_name=file[1]), CATEGORY, Qgis.Info)
                    reprojected_tiff = os.path.join(os.path.dirname(file[0]), "{filename}_NS_Corrected{ext}".format(filename=file[1],ext=os.path.splitext(file[0])[1]))
                    #Use gdalwarp to corrent the positive NS resolution
                    gdal.Warp(reprojected_tiff, tile_ds)
                    tile_ds = None

                    #Delete original file
                    os.remove(file[0])
                    #Rename reprojected tiff to original file name
                    os.rename(reprojected_tiff, file[0])
                else:
                    tile_ds = None

                org_files.append(file[0])
            
            except Exception as e:
                QgsMessageLog.logMessage('Skipping file: {1} | Error: {0}'.format(str(e), file[0]), CATEGORY, Qgis.Info)
                taskData.setTilesCount(0)
                return taskData
        
        org_vrt_options = None

        if manual_nodata_value is not None:
            org_vrt_options = gdal.BuildVRTOptions(resolution="highest",resampleAlg="bilinear",srcNodata=manual_nodata_value,VRTNodata=manual_nodata_value)
        else:
            nd_ds = gdal.Open(org_files[0], gdal.GA_ReadOnly)
            nd = nd_ds.GetRasterBand(1).GetNoDataValue()
            if nd is None:
                nd = 'nan'
            nd_ds = None
            
            org_vrt_options = gdal.BuildVRTOptions(resolution="highest",resampleAlg="bilinear", VRTNodata=nd)

        ds = gdal.BuildVRT(org_vrt, org_files,options=org_vrt_options)
        ds.FlushCache()

        source_no_data_value = None
        os_name = QgsApplication.osName()

        if os_name == "linux":
            if manual_nodata_value is None:
                huge_parser = etree.XMLParser(huge_tree=True)
                #Get no data value from org source vrt as string
                org_vrt_tree = etree.parse(org_vrt, parser=huge_parser)
                org_vrt_root = org_vrt_tree.getroot()
                org_vrt_band_1 = org_vrt_root.findall(".//VRTRasterBand[@band='1']")[0]
                org_vrt_band_1_no_data_entries = org_vrt_band_1.findall('.//NoDataValue')
                if org_vrt_band_1_no_data_entries is not None and len(org_vrt_band_1_no_data_entries) > 0:
                    org_vrt_band_1_no_data_entry = org_vrt_band_1_no_data_entries[0]
                    source_no_data_value = str(org_vrt_band_1_no_data_entry.text)

                if source_no_data_value == "" or source_no_data_value is None:
                    source_no_data_value = None
            else:
                source_no_data_value = str(manual_nodata_value)


        ds = None
        sleep(0.05)

        QgsMessageLog.logMessage(
                    'Created original vrt',
                    CATEGORY, Qgis.Info)

        vrt = os.path.join(source_folder, "Source.vrt").replace("\\","/")
        if os.path.isfile(vrt):
            os.remove(vrt)

        pds = gdal.Open(org_vrt)
        pds = gdal.Warp(vrt, pds, dstSRS="EPSG:3857")
        pds = None

        QgsMessageLog.logMessage('Created reprojected vrt',CATEGORY, Qgis.Info)

        terrarium_tile =  os.path.join(source_folder, "TerrariumSource.vrt").replace("\\","/")
        if os.path.isfile(terrarium_tile):
            os.remove(terrarium_tile)


        statistics = []

        task.setProgress(0)

        exxf = concurrent.futures.ThreadPoolExecutor(max_workers=cpu_count)
        fff= exxf.map(calculateStat, org_files)
        rt = 0
        cf = len(org_files)
        for res in fff:
            if res is not None:
                rt += 1
                statistics.append(res)
                task.setProgress(max(0, min(int((rt * 12) / cf), 100)))
            else:
                cf -= 1

        org_files = None

        # Find minV and maxV from statistics

        minV = None
        maxV = None

        if len(statistics) > 1:
            for stat in statistics:
                if stat.minV is not None and stat.maxV is not None:
                    if minV is None:
                        minV = stat.minV
                    elif stat.minV < minV:
                        minV = stat.minV

                    if maxV is None:
                        maxV = stat.maxV
                    elif stat.maxV > maxV:
                        maxV = stat.maxV
        else:
            minV = statistics[0].minV
            maxV = statistics[0].maxV

        if minV is None or maxV is None:
            QgsMessageLog.logMessage('Error: Minimum and maximum height are None',CATEGORY, Qgis.Info)
            return None
        else:
            QgsMessageLog.logMessage('Minimum and maximum height are {mn} and {ma}'.format(mn=minV,ma=maxV),CATEGORY, Qgis.Info)

        statistics = None

        color_ramp_file = os.path.join(source_folder, "TerrariumSource.txt").replace("\\","/")
        if os.path.isfile(color_ramp_file):
            os.remove(color_ramp_file)
        color_ramp = []
        altitudes = []

        # Create color ramp
        for meters in range(minV, maxV):
            for fraction in range (0, 256):
                altitudes.append(meters+(fraction/256))

        exx = concurrent.futures.ThreadPoolExecutor(max_workers=cpu_count)
        ff= exx.map(createRamp, repeat(convert_feet_to_meters), altitudes)
        for res in ff:
            if res is not None:
                color_ramp.append(res)


        QgsMessageLog.logMessage(
                    'Creating color ramp file',
                    CATEGORY, Qgis.Info)
        
        color_ramp.sort(key = lambda x: x.altitude)
        f = open(color_ramp_file, "w")
        rt = 0
        cf = len(color_ramp)
        for ramp in color_ramp:
            if manual_nodata_value is not None:
                if ramp.altitude == manual_nodata_value:
                    f.write('{altitude}\t0\t0\t0\t0\n'.format(altitude=manual_nodata_value))
                    continue
            f.write('{altitude}\t{red}\t{green}\t{blue}\n'.format(altitude=ramp.altitude, red=ramp.red, green=ramp.green, blue=ramp.blue))
            rt += 1
            task.setProgress(max(0, min(int(((rt * 8) / cf) + 12), 100)))
        f.write('nv\t0\t0\t0\t0')
        f.close()
        sleep(0.05)

        QgsMessageLog.logMessage(
                    'Created color ramp file',
                    CATEGORY, Qgis.Info)

        QgsMessageLog.logMessage(
                    'Rendering vrt to terrarium format. This might take a while',
                    CATEGORY, Qgis.Info)

        dst_ds = gdal.Open(vrt)
        ds = gdal.DEMProcessing(terrarium_tile, dst_ds, 'color-relief', colorFilename=color_ramp_file, format="VRT", addAlpha=True)
        sleep(0.05)
        dst_ds = None

        QgsMessageLog.logMessage(
                    'Created vrt in terrarium format',
                    CATEGORY, Qgis.Info)
        
        #gdal.DEMProcessing with the 'color-relief' option and when the output datasource is a VRT, 
        #inserts two key-pairs at the start of  of the LUT (Look up table) element of each VRTRasterBand (including the 4th, ergo Alpha band), which It's keys are very similar to the source NODATA value
        #It does this on both Windows and Linux, with the only difference being that on Windows these two keys are more similar,
        #if not the first key being the same as the source NODATA, but on Linux these two keys are less similar to the source NODATA, compared to Windows
        #This, alongside with difference between the first and second key being greater then on Windows, 
        #and the second value pointing to the same byte value as the 3rd key-pair value (The real min height value of the dataset),
        #makes gdal/QGIS render the NODATA value as half transparent (ex. 128 byte), making the script render tiles with half-transparent areas
        #
        #To fix this bug, the script first reads the key-pairs of the LUT element of the Alpha VRTRasterBand
        #If the first key-pair key is not the same as the NODATA value, it calculates the string similarity between the first and second key-pair keys
        #If it's the similarity ratio is greater then 0.9, the script removed the 2nd key-pair
        #Finally set's the first key-pair key to the NODATA value, It's value to 0 (0 byte, ergo Transparent) 
        #and writes the changes to the VRT file
        if os_name == "linux" and source_no_data_value is not None:
            QgsMessageLog.logMessage(
                    'Checking for no data value in alpha raster band of terrarium vrt',
                    CATEGORY, Qgis.Info)
            
            huge_parser = etree.XMLParser(huge_tree=True)
            terrarium_tile_tree = etree.parse(terrarium_tile,parser=huge_parser)
            terrarium_tile_root = terrarium_tile_tree.getroot()

            alpha_vrt_raster_band = terrarium_tile_root.findall(".//VRTRasterBand[@band='4']")[0]
            alpha_band_lut = alpha_vrt_raster_band.findall(".//LUT")[0]
            first_lut_key= alpha_band_lut.text[0:alpha_band_lut.text.index(',') - 2]
            if first_lut_key != source_no_data_value:
                
                org_lut_text = alpha_band_lut.text[alpha_band_lut.text.index(',') + 1:len(alpha_band_lut.text)]
                second_lut_key_index = org_lut_text.index(',')

                #Compare if the second LUT key is similar to the first one. If it is, remove it
                if second_lut_key_index > -1:
                    second_lut_key = org_lut_text[0:second_lut_key_index].split(':')[0]
                    first_second_match = SequenceMatcher(None, first_lut_key, second_lut_key)
                    if first_second_match.ratio() > 0.9:
                        org_lut_text = org_lut_text[second_lut_key_index + 1:len(org_lut_text)]
                
                alpha_band_lut_text = source_no_data_value + ':0,' + org_lut_text 
                alpha_band_lut.text = alpha_band_lut_text

                new_tree = etree.ElementTree(terrarium_tile_root)
                new_tree.write(terrarium_tile,pretty_print=True,xml_declaration=False,encoding="utf-8")

                QgsMessageLog.logMessage(
                    'Detected incorrect no data value ({0}) in alpha raster band of terrarium vrt. Should be: {1}. Corrected to right value'.format(first_lut_key, source_no_data_value),
                    CATEGORY, Qgis.Info)

        input_file = terrarium_tile

        QgsMessageLog.logMessage('Input file is {inp}'.format(inp=input_file), CATEGORY, Qgis.Info)
        dst_ds = gdal.Open(input_file, gdal.GA_ReadOnly)

        QgsMessageLog.logMessage('Calculating bounds', CATEGORY, Qgis.Info)

        diff = (20037508.3427892439067364 + 20037508.3427892439067364) / 2**zoom

        x_diff = diff
        y_diff = diff

        ulx, xres, xskew, uly, yskew, yres = dst_ds.GetGeoTransform()
        #lrx = ulx + (dst_ds.RasterXSize * xres)
        #lry = uly + (dst_ds.RasterYSize * yres)

        info = gdal.Info(dst_ds, format='json')

        ulx, uly = info['cornerCoordinates']['upperLeft'][0:2]
        lrx, lry = info['cornerCoordinates']['lowerRight'][0:2]

        wgs_minx, wgs_minz = info['wgs84Extent']['coordinates'][0][1][0:2]
        wgs_maxx, wgs_maxz = info['wgs84Extent']['coordinates'][0][3][0:2]

        if single_source_mode:
            QgsMessageLog.logMessage('The dataset bounds are (in WGS84 [EPSG:4326]), minX: {minX}, minZ: {minZ}, maxX: {maxX}, maxZ: {maxZ}'.format(minX=str(wgs_minx),minZ=str(wgs_minz),maxX=str(wgs_maxx),maxZ=str(wgs_maxz)), CATEGORY, Qgis.Info)
        else:
            QgsMessageLog.logMessage('The dataset bounds for the {dataset_name} dataset are (in WGS84 [EPSG:4326]), minX: {minX}, minZ: {minZ}, maxX: {maxX}, maxZ: {maxZ}'.format(dataset_name=taskData.getDatasetName(), minX=str(wgs_minx),minZ=str(wgs_minz),maxX=str(wgs_maxx),maxZ=str(wgs_maxz)), CATEGORY, Qgis.Info)
        
        QgsMessageLog.logMessage('Optional: Take note of these bounds if you want to manually create the entries in heights.json \n as described in step F. of Part two: Generating/using your dataset',CATEGORY, Qgis.Info)

        taskData.setSourceBounds(wgs_minx, wgs_minz, wgs_maxx, wgs_maxz)
        
        base_x = 0
        base_y = 0

        start_tile_x = base_x
        start_tile_y = base_y
        min_x = -20037508.3427892439067364
        max_y = 20037508.3427892439067364

        x_tiles = 0
        y_tiles = 0

        # Find min lot
        lon = min_x
        while lon <= 20037508.3427892439067364:
            if lon >= ulx:
                break;
            start_tile_x += 1
            min_x = lon
            lon += x_diff

        # Find max lat
        lat = max_y
        while lat >= -20037508.3427892550826073:
            if lat <= uly:
                break
            start_tile_y += 1
            max_y = lat
            lat -= y_diff

        # Find how many lon tiles to make
        lon = min_x
        while lon < lrx:
            x_tiles += 1
            lon += x_diff

        # Find how many lat tiles to make
        lat = max_y
        while lat >= lry:
            y_tiles += 1
            lat -= y_diff

        if start_tile_x > 0:
            start_tile_x -= 1
        if start_tile_y > 0:
            start_tile_y -= 1

        QgsMessageLog.logMessage('Start tile: {tx} ({mlon}), {ty} ({mxlat}), Tiles to generate: {xt} (Width: {xtx} Height: {xty})'.format(tx=start_tile_x, mlon=min_x, ty=start_tile_y, mxlat=max_y, xt=((x_tiles + 1) * (y_tiles + 1)), xtx=x_tiles,xty=y_tiles), CATEGORY, Qgis.Info)

        if single_source_mode:
            QgsMessageLog.logMessage('Creating output folders', CATEGORY, Qgis.Info)
        else:
            QgsMessageLog.logMessage('Creating output folders for {dataset_name} dataset'.format(dataset_name=taskData.getDatasetName()), CATEGORY, Qgis.Info)

        mx = start_tile_x + x_tiles
        zoom_folder = ''
        temp_folder = os.path.join(source_folder, 'TEMP')
        if not ftp_upload:
            if single_source_mode:
                zoom_folder = os.path.join(output, str(zoom)).replace("\\","/")
            else:
                output_dataset = os.path.join(output, taskData.getDatasetName()).replace("\\","/")
                if not os.path.isdir(output_dataset):
                    os.mkdir(output_dataset)

                zoom_folder = os.path.join(output_dataset, str(zoom)).replace("\\","/")
        else:
            if os.path.isdir(temp_folder):
                shutil.rmtree(temp_folder)
            os.mkdir(temp_folder)
            zoom_folder = os.path.join(temp_folder ,str(zoom)).replace("\\","/")
        if not os.path.isdir(zoom_folder):
            os.mkdir(zoom_folder)
        for x in range(start_tile_x, mx):
            folderx = os.path.join(zoom_folder, str(x)).replace("\\","/")
            if not os.path.isdir(folderx):
                os.mkdir(folderx)

        zip_file_name = 'RenderedDataset.zip'
        if single_source_mode is False:
            zip_file_name = 'RenderedDataset_{dataset_name}.zip'.format(dataset_name=taskData.getDatasetName())
        
        zip_file = os.path.join(zoom_folder, zip_file_name).replace("\\","/")
        rd_file = None
        if ftp_upload:
            if ftp_one_file:
                rd_file = zipfile.ZipFile(zip_file, 'w', zipfile.ZIP_DEFLATED)
            else:
                QgsMessageLog.logMessage('Creating output folders on target ftp server', CATEGORY, Qgis.Info)
                ftp = None
                if ftp_s:
                    ftp = FTP_TLS()
                else:
                    ftp = FTP()
                ftp.connect(ftp_upload_url, ftp_upload_port)
                if ftp_user is None or ftp_password is None:
                    ftp.login()
                else:
                    ftp.login(user=ftp_user, passwd=ftp_password)

                if output is not None:
                    ftp.cwd(output)
                
                if single_source_mode is False:
                    if not isFTPDir(ftp, taskData.getDatasetName()):
                        ftp.mkd(taskData.getDatasetName())
                    ftp.cwd(taskData.getDatasetName())


                if not isFTPDir(ftp, str(zoom)):
                    ftp.mkd(str(zoom))

                ftp.cwd(str(zoom))

                for x in range(start_tile_x, mx):
                    ftp.mkd(str(x))

                ftp.quit()


        sleep(0.01)

        QgsMessageLog.logMessage('Created {ct} folders'.format(ct=x_tiles), CATEGORY, Qgis.Info)

        tiled = []

        # Tile dataset
        sub_min_x = min_x
        sub_max_x = min_x + x_diff
        for x in range(start_tile_x, start_tile_x + x_tiles):
            sub_min_y = max_y - y_diff
            sub_max_y = max_y
            for y in range(start_tile_y, start_tile_y + y_tiles):
                tiled.append(Tile(dst_ds, str(x), str(y), sub_min_x, sub_max_y, sub_max_x, sub_min_y, ulx, xres, uly, yres, querysize=1024))
                sub_min_y -= y_diff
                sub_max_y -= y_diff
            sub_min_x += x_diff
            sub_max_x += x_diff

        job = Job(zoom_folder, str(zoom), input_file, getBands(dst_ds), resampling_algorithm, output, ftp_upload, ftp_one_file, ftp_s ,ftp_upload_url, ftp_upload_port, ftp_user, ftp_password, taskData.getDatasetName())

        dst_ds = None

        # Tile the dataset

        realtiles = 0

        if cpu_count == 1:
            if single_source_mode:
                QgsMessageLog.logMessage('Started tilling vrt in single-thread mode', CATEGORY, Qgis.Info)
            else:
                QgsMessageLog.logMessage('Started tilling vrt for {dataset_name}({index}/{length}) dataset in single-thread mode'.format(dataset_name=taskData.getDatasetName(), index=taskData.sourceIndex + 1, length=sources_length), CATEGORY, Qgis.Info)
            rt = 0
            cf = len(tiled)
            for tile in tiled:
                res = tileVrt(job, tile)
                if res is not None:
                    if job.ftpUpload and job.ftpOnefile:
                        addToZIP(job, res, rd_file)
                    realtiles += 1
                rt += 1
                if job.ftpUpload and job.ftpOnefile:
                    task.setProgress(max(0, min(int(((rt * 40) / cf) + 20), 100)))
                else:
                    task.setProgress(max(0, min(int(((rt * 80) / cf) + 20), 100)))

            if getattr(localThread, 'ds', None):
                del localThread.ds

        else:
            if single_source_mode:
                QgsMessageLog.logMessage('Started tilling vrt in multithread mode with {count} threads'.format(count=cpu_count), CATEGORY, Qgis.Info)
            else:
                QgsMessageLog.logMessage('Started tilling vrt for {dataset_name}({index}/{length}) dataset in multithread mode with {count} threads'
                                         .format(dataset_name=taskData.getDatasetName(), index=taskData.sourceIndex + 1, length=sources_length,count=cpu_count), CATEGORY, Qgis.Info)

            gdal_cache_max = gdal.GetCacheMax()
            gdal_cache_max_per_process = max(1024 * 1024, math.floor(gdal_cache_max / cpu_count))
            setCacheMax(gdal_cache_max_per_process)


            tpe = concurrent.futures.ThreadPoolExecutor(max_workers=cpu_count)
            tm = tpe.map(tileVrt, repeat(job), tiled)
            rt = 0
            cf = len(tiled)
            realtiles = 0
            for res in tm:
                if res is not None:
                    if job.ftpUpload and job.ftpOnefile:
                        addToZIP(job, res, rd_file)
                    realtiles += 1
                rt += 1
                if job.ftpUpload and job.ftpOnefile:
                    task.setProgress(max(0, min(int(((rt * 40) / cf) + 20), 100)))
                else:
                    task.setProgress(max(0, min(int(((rt * 80) / cf) + 20), 100)))

            setCacheMax(gdal_cache_max)

        if job.ftpUpload and job.ftpOnefile:
            rd_file.close()
            totalSize = os.path.getsize(zip_file)
            QgsMessageLog.logMessage('Starting uploading rendered archive ({size} GB) out of {count} tiles'.format(count=realtiles, size=str(round(totalSize * 9.3132257461548E-10, 4))), CATEGORY, Qgis.Info)

            ftp = None
            if job.ftpS:
                ftp = FTP_TLS()
            else:
                ftp = FTP()
            ftp.connect(job.ftpUrl, job.ftpPort)

            if job.ftpUser is None or job.ftpPassword is None:
                ftp.login()
            else:
                ftp.login(user=job.ftpUser, passwd=job.ftpPassword)

            if job.outputPath is not None:
                ftp.cwd(job.outputPath)

            ftapr = FTPArchiveProgress(int(totalSize), task)

            with open(zip_file, 'rb') as rdf_file:
                ftp.storbinary('STOR {dst_zip_file}'.format(dst_zip_file=zip_file_name), rdf_file, blocksize=1024, callback=ftapr.updateProgress)

            ftp.quit()
            
            if single_source_mode:
                QgsMessageLog.logMessage('Uploaded rendered dataset archive {dst_zip_file} with {count} tiles. You can now unzip it.'.format(dst_zip_file=zip_file_name, count=realtiles), CATEGORY, Qgis.Info)
            else:
                QgsMessageLog.logMessage('Uploaded rendered dataset archive {dst_zip_file}({index}/{length}) with {count} tiles. You can now unzip it.'.format(dst_zip_file=zip_file_name,index=taskData.sourceIndex + 1, length=sources_length, count=realtiles), CATEGORY, Qgis.Info)
        elif job.ftpUpload:
            if single_source_mode:
                QgsMessageLog.logMessage('Tiled and uploaded dataset with {count} tiles to ftp server'.format(count=realtiles), CATEGORY, Qgis.Info)
            else:
                QgsMessageLog.logMessage('Tiled and uploaded {dataset_name} dataset ({index}/{length})  with {count} tiles to ftp server'.format(dataset_name=taskData.getDatasetName(), index=taskData.sourceIndex + 1, length=sources_length, count=realtiles), CATEGORY, Qgis.Info)
        else:
            if single_source_mode:
                QgsMessageLog.logMessage('Tiled dataset with {count} tiles'.format(count=realtiles), CATEGORY, Qgis.Info)
            else:
                QgsMessageLog.logMessage('Tiled {dataset_name} dataset ({index}/{length}) with {count} tiles'.format(dataset_name=taskData.getDatasetName(), index=taskData.sourceIndex + 1, length=sources_length, count=realtiles), CATEGORY, Qgis.Info)

        # Clean up
        if cleanup:
            os.remove(org_vrt)
            os.remove(vrt)
            os.remove(color_ramp_file)
            os.remove(terrarium_tile)

            if job.ftpUpload:
                shutil.rmtree(temp_folder)

            QgsMessageLog.logMessage('Cleaned up temp files', CATEGORY, Qgis.Info)


        taskData.setTilesCount(realtiles)
        return taskData
    except Exception as e:
        QgsMessageLog.logMessage('Error: ' + str(e), CATEGORY, Qgis.Info)
        traceback.print_exc()
        return None


def setCacheMax(cache_in_bytes: int):
    os.environ['GDAL_CACHEMAX'] = '%d' % int(cache_in_bytes / 1024 / 1024)
    gdal.SetCacheMax(cache_in_bytes)

def getBands(src_ds):
    bandsCount = 0
    alphaband = src_ds.GetRasterBand(1).GetMaskBand()
    if ((alphaband.GetMaskFlags() & gdal.GMF_ALPHA) or
            src_ds.RasterCount == 4 or
            src_ds.RasterCount == 2):
        bandsCount = src_ds.RasterCount - 1
    else:
        bandsCount = src_ds.RasterCount
    return bandsCount

def createRamp(cftm, altitude):
        return ColorRamp(altitude, cftm)

def calculateStat(tile):
    try:
        stat_ds = gdal.Open(tile)

        band = stat_ds.GetRasterBand(1)

        minV = None
        maxV = None

        if band.GetMinimum() is None or band.GetMaximum() is None:
            stats = band.GetStatistics(0,1)
            if band.GetMinimum() is not None and band.GetMaximum() is not None:
                minV = int(math.floor(stats[0]))
                maxV = int(math.ceil(stats[1]))
            else:
                stat_ds = None
                QgsMessageLog.logMessage(
                            'Tile: {name} is empty. Ignoring stats.'.format(name=tile),
                            CATEGORY, Qgis.Info)   
                return None
        else:
            minV = int(math.floor(band.GetMinimum()))
            maxV = int(math.ceil(band.GetMaximum()))

        stat_ds = None

        return Statistic(minV, maxV)
    except Exception as e:
        QgsMessageLog.logMessage(
                    'Error while calculating stat for tile: {name}. Error: {er}'.format(name=tile,er=str(e)),
                    CATEGORY, Qgis.Info)
        return None

def isFTPDir(ftp, name):
   try:
      ftp.cwd(name)
      ftp.cwd('..')
      return True
   except:
      return False

def addToZIP(job_data, tile_data, rd_file):
    folder = os.path.join(job_data.zoom_folder, tile_data.x).replace("\\","/")
    tileName = tile_data.y + ".png"
    tile_file =  os.path.join(folder, tileName).replace("\\","/")

    #Add tile to zip
    rd_file.write(tile_file, os.path.join(job_data.zoom_level, str(tile_data.x), str(tile_data.y) + ".png"))
    #Remove tile from TEMP
    os.remove(tile_file)

# Modified version of scale_query_to_tile from gdal2tiles
#https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/scripts/gdal2tiles.py
def scaleTile(dsquery, dstile, resampling, tile=''):
    querysize = dsquery.RasterXSize
    tile_size = dstile.RasterXSize
    tilebands = dstile.RasterCount

    if resampling == 'average':

        # Function: gdal.RegenerateOverview()
        for i in range(1, tilebands + 1):
            # Black border around NODATA
            res = gdal.RegenerateOverview(dsquery.GetRasterBand(i), dstile.GetRasterBand(i),
                                          'average')
            if res != 0:
                QgsMessageLog.logMessage("RegenerateOverview() failed on %s, error %d" % (
                    tile, res), CATEGORY, Qgis.Info)

    elif resampling == 'antialias' and numpy_available:

        # Scaling by PIL (Python Imaging Library) - improved Lanczos
        array = numpy.zeros((querysize, querysize, tilebands), numpy.uint8)
        for i in range(tilebands):
            array[:, :, i] = gdalarray.BandReadAsArray(dsquery.GetRasterBand(i + 1),
                                                       0, 0, querysize, querysize)
        im = Image.fromarray(array, 'RGBA')     # Always four bands
        im1 = im.resize((tile_size, tile_size), Image.ANTIALIAS)
        if os.path.exists(tile):
            im0 = Image.open(tile)
            im1 = Image.composite(im1, im0, im1)
        im1.save(tile, 'PNG')

    else:

        if resampling == 'near':
            gdal_resampling = gdal.GRA_NearestNeighbour

        elif resampling == 'bilinear':
            gdal_resampling = gdal.GRA_Bilinear

        elif resampling == 'cubic':
            gdal_resampling = gdal.GRA_Cubic

        elif resampling == 'cubicspline':
            gdal_resampling = gdal.GRA_CubicSpline

        elif resampling == 'lanczos':
            gdal_resampling = gdal.GRA_Lanczos

        elif resampling == 'mode':
            gdal_resampling = gdal.GRA_Mode

        elif resampling == 'max':
            gdal_resampling = gdal.GRA_Max

        elif resampling == 'min':
            gdal_resampling = gdal.GRA_Min

        elif resampling == 'med':
            gdal_resampling = gdal.GRA_Med

        elif resampling == 'q1':
            gdal_resampling = gdal.GRA_Q1

        elif resampling == 'q3':
            gdal_resampling = gdal.GRA_Q3

        # Other algorithms are implemented by gdal.ReprojectImage().
        dsquery.SetGeoTransform((0.0, tile_size / float(querysize), 0.0, 0.0, 0.0,
                                 tile_size / float(querysize)))
        dstile.SetGeoTransform((0.0, 1.0, 0.0, 0.0, 0.0, 1.0))

        res = gdal.ReprojectImage(dsquery, dstile, None, None, gdal_resampling)
        if res != 0:
            QgsMessageLog.logMessage("ReprojectImage() failed on %s, error %d" % (tile, res), CATEGORY, Qgis.Info)


# Modified version of create_base_tile from gdal2tiles
#https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/scripts/gdal2tiles.py
def tileVrt(job_data, tile_data):
    try:
        folder = os.path.join(job_data.zoom_folder, tile_data.x).replace("\\","/")
        tileName = tile_data.y + ".png"
        tile =  os.path.join(folder, tileName).replace("\\","/")

        tilebands = job_data.bandsCount + 1

        ds_cache = getattr(localThread, 'ds', None)
        if ds_cache:
            ds = ds_cache
        else:
            ds = gdal.Open(job_data.input_file, gdal.GA_ReadOnly)
            localThread.ds = ds

        mem_drv = gdal.GetDriverByName('MEM')
        out_drv = gdal.GetDriverByName('PNG')
        alphaband = ds.GetRasterBand(1).GetMaskBand()

        dstile = mem_drv.Create('', 256, 256, job_data.bandsCount + 1)
        data = alphamask = None

        if tile_data.rxsize != 0 and tile_data.rysize != 0 and tile_data.wxsize != 0 and tile_data.wysize != 0:
            alphamask = alphaband.ReadRaster(tile_data.rx, tile_data.ry, tile_data.rxsize, tile_data.rysize, tile_data.wxsize, tile_data.wysize)

            ##Check if transparent
            if len(alphamask) == alphamask.count('\x00'.encode('ascii')):
                return None

            data = ds.ReadRaster(tile_data.rx, tile_data.ry, tile_data.rxsize, tile_data.rysize, tile_data.wxsize, tile_data.wysize, band_list=list(range(1, job_data.bandsCount + 1)))

        if data:
            if  tile_data.querysize == 256:
                dstile.WriteRaster(tile_data.wx, tile_data.wy, tile_data.wxsize, tile_data.wysize, data,
                                   band_list=list(range(1, job_data.bandsCount + 1)))
                dstile.WriteRaster(tile_data.wx, tile_data.wy, tile_data.wxsize, tile_data.wysize, alphamask, band_list=[tilebands])
            else:
                dsquery = mem_drv.Create('', tile_data.querysize, tile_data.querysize, tilebands)
                dsquery.WriteRaster(tile_data.wx, tile_data.wy, tile_data.wxsize, tile_data.wysize, data,
                                    band_list=list(range(1, job_data.bandsCount + 1)))
                dsquery.WriteRaster(tile_data.wx, tile_data.wy, tile_data.wxsize, tile_data.wysize, alphamask, band_list=[tilebands])

                scaleTile(dsquery, dstile, job_data.resampling, tile=tile)
                del dsquery

        del data

        if job_data.resampling != 'antialias':
            out_drv.CreateCopy(tile, dstile, strict=0)

        del dstile

        xmlFile = os.path.join(folder, tile_data.y + '.png.aux.xml')
        if os.path.isfile(xmlFile):
            os.remove(xmlFile)


        if job_data.ftpUpload and not job_data.ftpOnefile:
            ftp = None
            if job_data.ftpS:
                ftp = FTP_TLS()
            else:
                ftp = FTP()
            ftp.connect(job_data.ftpUrl, job_data.ftpPort)

            if job_data.ftpUser is None or job_data.ftpPassword is None:
                ftp.login()
            else:
                ftp.login(user=job_data.ftpUser, passwd=job_data.ftpPassword)

            if job_data.outputPath is not None:
                ftp.cwd(job_data.outputPath)

            if single_source_mode is False:
                ftp.cwd(job_data.getDatasetName())
            
            ftp.cwd(str(job_data.zoom_level) + '/' + tile_data.x)
            with open(tile, 'rb') as tile_file:
                ftp.storbinary('STOR {tileName}'.format(tileName=tileName), tile_file)

            ftp.quit()
            os.remove(tile)


        sleep(0.01)
    except Exception as e:
        QgsMessageLog.logMessage('Tiled vrt error: ' + str(e), CATEGORY, Qgis.Info)
        traceback.print_exc()
        return None

    return tile_data

def getHeightsEntry(url_path, zoom_level, minx, minz, maxx, maxz):
    height_entry = {
        "dataset" : {
            "urls" : [ url_path],
            "projection" : {
                "web_mercator" : {
                    "zoom" : zoom_level
                }
            },
            "resolution" : 256,
            "blend" : "CUBIC",
            "parse" : {
                "parse_png_terrarium" : {}
            }
        },
        "bounds" : {
            "minX" : minx,
            "maxX" : maxx,
            "minZ" : minz,
            "maxZ" : maxz
        },
        "zooms" : {
            "min" : 0,
            "max" : 3
        },
        "priority" : 100
    }

    return height_entry
    

def genHeightsConfig(taskData):
    sources_array = taskData.sourcesArray
    heights_list = []

    for source in sources_array:
        wgs_minx, wgs_minz, wgs_maxx, wgs_maxz = source.sourceBounds
        url_path = "file://"
        if output is not None:
            url_path += f"{output}/"
        if single_source_mode is False:
            url_path += f"{source.datasetName}/"
        url_path += str(zoom) + "/${x}/${z}.png"

        heights_list.append(getHeightsEntry(url_path, zoom, wgs_minx, wgs_minz, wgs_maxx, wgs_maxz))

    raw_json = json.dumps(heights_list, indent=4)

    if ftp_upload:
        #Upload heights config to the server via ftp
        
        ftp = None
        if ftp_s:
            ftp = FTP_TLS()
        else:
            ftp = FTP()
        ftp.connect(ftp_upload_url, ftp_upload_port)

        if ftp_user is None or ftp_password is None:
            ftp.login()
        else:
            ftp.login(user=ftp_user, passwd=ftp_password)

        if output is not None:
            ftp.cwd(output)

        json_byte = io.BytesIO()
        json_byte.write(raw_json.encode())
        json_byte.seek(0)

        ftp.storbinary('STOR heightsTemplate.json', json_byte)

        if output is not None:
            QgsMessageLog.logMessage(f"Uploaded heightsTemplate.json to the {output} folder on your server", 
                CATEGORY, 
                Qgis.Info)
        else:
            QgsMessageLog.logMessage("Uploaded heightsTemplate.json to the root folder of your server", 
                CATEGORY, 
                Qgis.Info)

    elif output is not None:
        heights_template_path = os.path.join(output, "heightsTemplate.json")

        with open(heights_template_path, "w") as heightsTemplate:
            heightsTemplate.write(raw_json)

        QgsMessageLog.logMessage(f"Created heightsTemplate.json in the {output} folder", 
                CATEGORY, 
                Qgis.Info)
    
    QgsMessageLog.logMessage("Head to step F. of Part two: Generating/using your dataset to set it up", CATEGORY, Qgis.Info)


def taskComplete(task, taskData=None):
    if taskData is not None:
        time_end = datetime.now()
        eclipsed = (time_end - taskData.getTimeStarted()).total_seconds() / 60.0
        minutes = math.floor(eclipsed)
        seconds = math.floor((eclipsed - minutes) * 60)
        if single_source_mode:
            QgsMessageLog.logMessage('Done creating dataset with {count} tiles in {minutes} minutes and {seconds} seconds'.format(count=taskData.getTilesCount(), minutes=minutes, seconds=seconds), CATEGORY, Qgis.Info)
        else:
            QgsMessageLog.logMessage('Done creating {dataset_name} dataset ({index}/{lenght}) with {count} tiles in {minutes} minutes and {seconds} seconds'.format(dataset_name=taskData.getDatasetName(), index=taskData.sourceIndex + 1, lenght=len(taskData.sourcesArray), count=taskData.getTilesCount(), minutes=minutes, seconds=seconds), CATEGORY, Qgis.Info)

            #Move onto next source in the sources array
            if taskData.setSourceIndex(taskData.sourceIndex + 1):
                globals()['dataset_task'] = QgsTask.fromFunction(taskData.genTaskName(), genTiles, on_finished=taskComplete, taskData=taskData)
                QgsApplication.taskManager().addTask(globals()['dataset_task'])                
            else:
                time_end = datetime.now()
                eclipsed = (time_end - taskData.timeStartedFirst).total_seconds() / 60.0
                minutes = math.floor(eclipsed)
                seconds = math.floor((eclipsed - minutes) * 60)
                QgsMessageLog.logMessage('Done creating all datasets in {minutes} minutes and {seconds} seconds'.format(minutes=minutes, seconds=seconds), CATEGORY, Qgis.Info)
                genHeightsConfig(taskData)

sources_list = []

for source in sources:
    source_path = source.replace("\\","/")
    if source_path not in sources_list:
        sources_list.append(source_path)

single_source_mode = True
if len(sources_list) > 1:
    single_source_mode = False

sources_array = []
for source in sources_list:
    sources_array.append(SourceDataset(source))

source_one_name = sources_array[0].datasetName
task_data = TaskData(0, sources_array, datetime.now())

task_one_name = 'Create elevation dataset'
if single_source_mode is False:
    task_one_name = task_data.genTaskName()

output = output.replace("\\","/")

globals()['dataset_task'] = QgsTask.fromFunction(task_one_name, genTiles, on_finished=taskComplete, taskData=task_data)
QgsApplication.taskManager().addTask(globals()['dataset_task'])