gdal_retile.py translated into C# beta

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Drawing;
using System.IO;
using System.Linq;
using System.Windows.Media.Imaging;
using System.Text;
using OSGeo.GDAL;
using OSGeo.OGR;
using OSGeo.OSR;
using Driver = OSGeo.GDAL.Driver;

namespace RasterLoader
{
    public class GdalRetile
    {
        private const double MetersPerInch = 0.0254;
        private const double DPI = 96;

        private string _format = string.Empty;
        private List<string> _createOptions = new List<string>();
        private DataType _bandType;
        private bool _verbose;
        private string _targetDir = string.Empty;
        private int _tileWidth = -1;
        private int _tileHeight = -1;
        private string _resamplingMethodString = string.Empty;
        private ResampleAlg _resamplingMethod;
        private int _levels = -1;
        private SpatialReference _sourceSRS;
        private bool _pyramidOnly;
        private string _tileIndexName = string.Empty;
        private string _tileIndexFieldName = string.Empty;
        private string _tileIndexDriverTyp = string.Empty;
        private string _csvFileName = string.Empty;
        private string _csvDelimiter = string.Empty;
        private bool _useDirForEachRow;
        private List<string> _names = new List<string>();
        private Driver _driver;
        private Driver _memDriver;
        private int _lastRowIndex;
        private string _extension;
        private bool _8to32bit;
        private double _scaleFactor;
        private double _rasterSize = -1;
        private double _meters = -1;

        public GdalRetile(string[] args)
        {
            Gdal.AllRegister();

            InitFields();

            if (args.Any())
            {
                string[] argv = Gdal.GeneralCmdLineProcessor(args, 0);
                if (!argv.Any())
                    throw new Exception("Command line has no parameters.");

                int i = 1;
                while (i < argv.Count())
                {
                    string arg = argv.ElementAt(i);
                    string ext;
                    switch (arg)
                    {
                        case "-of":
                            i++;
                            _format = argv.ElementAt(i);
                            break;
                        case "-ot":
                            i++;
                            _bandType = Gdal.GetDataTypeByName(argv.ElementAt(i));
                            Console.WriteLine(_bandType);
                            if (_bandType == DataType.GDT_Unknown)
                            {
                                throw new Exception("Unknown GDAL data type: " + argv.ElementAt(i));
                            }
                            break;
                        case "-co":
                            i++;
                            _createOptions.Add(argv.ElementAt(i));
                            break;
                        case "-v":
                            _verbose = true;
                            break;
                        case "-targetDir":
                            i++;
                            _targetDir = argv.ElementAt(i);
                            if (!Directory.Exists(_targetDir))
                            {
                                throw new Exception("TargetDir " + _targetDir + " does not exist.");
                            }
                            break;
                        case "-ps":
                            i++;
                            _tileWidth = int.Parse(argv.ElementAt(i));
                            i++;
                            _tileHeight = int.Parse(argv.ElementAt(i));
                            break;
                        case "-r":
                            i++;
                            _resamplingMethodString = argv.ElementAt(i);
                            switch (_resamplingMethodString)
                            {
                                case "near":
                                    _resamplingMethod = ResampleAlg.GRA_NearestNeighbour;
                                    break;
                                case "bilinear":
                                    _resamplingMethod = ResampleAlg.GRA_Bilinear;
                                    break;
                                case "cubic":
                                    _resamplingMethod = ResampleAlg.GRA_Cubic;
                                    break;
                                case "cubicspline":
                                    _resamplingMethod = ResampleAlg.GRA_CubicSpline;
                                    break;
                                case "lanczos":
                                    _resamplingMethod = (ResampleAlg)4;
                                    break;
                                default:
                                    throw new Exception("Unknown resampling method: " + _resamplingMethodString);
                            }
                            break;
                        case "-levels":
                            i++;
                            _levels = int.Parse(argv.ElementAt(i));
                            if (_levels < 1)
                            {
                                throw new Exception("Invalid number of _levels: " + _levels);
                            }
                            break;
                        case "-s_srs":
                            i++;
                            string wkt = string.Empty;
                            if (Osr.GetUserInputAsWKT(argv.ElementAt(i), out wkt) != 0)
                            {
                                throw new Exception("Invalid -s_srs: " + argv.ElementAt(i));
                            }
                            _sourceSRS = new SpatialReference(wkt); // TODO: check this!
                            break;
                        case "-pyramidOnly":
                            _pyramidOnly = true;
                            break;
                        case "-tileIndex":
                            i++;
                            _tileIndexName = argv.ElementAt(i);
                            ext = Path.GetExtension(_tileIndexName);
                            if (ext.Length == 0)
                            {
                                _tileIndexName += ".shp";
                            }
                            break;
                        case "-tileIndexField":
                            i++;
                            _tileIndexFieldName = argv.ElementAt(i);
                            break;
                        case "-csv":
                            i++;
                            _csvFileName = argv.ElementAt(i);
                            ext = Path.GetExtension(_csvFileName);
                            if (ext.Length == 0)
                            {
                                _csvFileName += ".csv";
                            }
                            break;
                        case "-csvDelim":
                            i++;
                            _csvDelimiter = argv.ElementAt(i);
                            break;
                        case "-useDirForEachRow":
                            _useDirForEachRow = true;
                            break;
                        case "-8to32bit":
                            _8to32bit = true;
                            break;
                        default:
                            if (arg.StartsWith("-"))
                            {
                                Usage();
                                throw new Exception("Unrecognized command option: " + arg);
                            }
                            _names.Add(arg);
                            break;
                    }
                    i++;
                }

                if (!_names.Any())
                {
                    Usage();
                    throw new Exception("No input files selected.");
                }

                if (_tileWidth <= 0 || _tileHeight <= 0)
                {
                    Usage();
                    throw new Exception("Invalid tile dimension: " + _tileWidth + "x" + _tileHeight + ".");
                }

                if (string.IsNullOrEmpty(_targetDir) || string.IsNullOrWhiteSpace(_targetDir))
                {
                    Usage();
                    throw new Exception("Missing directory for tiles [-targetDir].");
                }

                if (_8to32bit)
                {
                    foreach (string name in _names)
                    {
                        //TiffConverter.Convert8To32Bit(name);
                    }
                }

                // create level 0 directory if needed
                if (_useDirForEachRow && _pyramidOnly == false)
                {
                    string levelDir = Path.Combine(_targetDir, 0.ToString());
                    if (!Directory.Exists(levelDir))
                        Directory.CreateDirectory(levelDir);
                }

                // prepare dirs for the pyramid
                if (_levels > 0)
                {
                    int startIndex = 1;
                    for (int levelIndex = startIndex; levelIndex < _levels + 1; levelIndex++)
                    {
                        string levelDir = Path.Combine(_targetDir, levelIndex.ToString());
                        if (Directory.Exists(levelDir))
                        {
                            continue;
                        }
                        try
                        {
                            Directory.CreateDirectory(levelDir);
                        }
                        catch
                        {
                            throw new Exception("Cannot create level directory: " + levelDir);
                        }
                        if (_verbose)
                        {
                            Debug.WriteLine("Created level directory: " + levelDir);
                        }
                    }
                }

                _driver = Gdal.GetDriverByName(_format);
                if (_driver == null)
                {
                    UsageFormat();
                    throw new Exception("Format _driver " + _format + " not found, pick a supported _driver.");
                }

                //_extension = _driver.GetMetadataItem(GdalConst.GDAL_DMD_EXTENSION, GdalConst.GDAL_DMD_LONGNAME); // TODO: check this
                string[] driverMD = _driver.GetMetadata(null); // null for the default domain
                _extension = _driver.GetMetadataItem("DMD_EXTENSION", null);

                if (driverMD.Contains("DCAP_CREATE"))
                {
                    _memDriver = Gdal.GetDriverByName("MEM");
                }

                DataSource tileIndexDS = GetTileIndexFromFiles(_names, _tileIndexDriverTyp);
                if (tileIndexDS == null)
                {
                    throw new Exception("Error building tile index.");
                }

                MosaicInfo mInfo = new MosaicInfo(_names[0], tileIndexDS);
                TileInfo ti = new TileInfo(mInfo.XSize, mInfo.YSize, _tileWidth, _tileHeight);

                if (_sourceSRS == null && mInfo.Projection.Length > 0)
                {
                    // TODO: check this
                    string wkt = string.Empty;
                    Osr.GetUserInputAsWKT(mInfo.Projection, out wkt);
                    _sourceSRS = new SpatialReference(wkt);

                    if (_sourceSRS.SetFromUserInput(mInfo.Projection) != 0)
                    {
                        throw new Exception("Invalid projection: " + mInfo.Projection);
                    }
                }

                if (_verbose)
                {
                    ti.Report();
                }

                DataSource dsCreatedTileIndex;
                if (!_pyramidOnly)
                {
                    dsCreatedTileIndex = TileImage(mInfo, ti);
                    tileIndexDS.Dispose();
                }
                else
                {
                    dsCreatedTileIndex = tileIndexDS;
                }

                if (_levels > 0)
                {
                    BuildPyramid(mInfo, dsCreatedTileIndex, _tileWidth, _tileHeight);
                }
            }
            else
            {
                throw new Exception("Command line has no parameters.");
            }
        }

        private void InitFields()
        {
            _verbose = false;
            _createOptions = new List<string>();
            _names = new List<string>();
            _tileWidth = 256;
            _tileHeight = 256;
            _format = "GTiff";
            _bandType = DataType.GDT_Unknown;
            _driver = null;
            _extension = null;
            _memDriver = null;
            _tileIndexFieldName = "location";
            _tileIndexName = null;
            _tileIndexDriverTyp = "Memory";
            _csvDelimiter = ";";
            _csvFileName = null;
            _sourceSRS = null;
            _targetDir = null;
            _resamplingMethod = ResampleAlg.GRA_NearestNeighbour;
            _levels = 0;
            _pyramidOnly = false;
            _lastRowIndex = -1;
            _useDirForEachRow = false;
        }

        private double GetScale(int meters, int pixels)
        {
            return meters / pixels / MetersPerInch * DPI;
        }

        private void BuildPyramid(MosaicInfo mInfo, DataSource createdTileIndexDS, int tileWidth, int tileHeight)
        {
            DataSource inputDS = createdTileIndexDS;

            for (int level = 1; level < _levels + 1; level++)
            {
                _lastRowIndex = -1;
                MosaicInfo levelMosaicInfo = new MosaicInfo(mInfo.FileName, inputDS);

                if (level == 1)
                {
                    double scale = GetScale((int) Math.Round(_meters), (int) Math.Round(_rasterSize));
                    double power = Math.Log(scale, 2);
                    if (power - (double) ((int) power) > 0d)
                    {
                        power = Math.Ceiling(power);
                    }

                    _scaleFactor = Math.Pow(2, power) / scale;
                }
                else
                {
                    _scaleFactor = 2;
                }
                TileInfo levelOutputTileInfo = new TileInfo(levelMosaicInfo.XSize / _scaleFactor, levelMosaicInfo.YSize / _scaleFactor, tileWidth, tileHeight);
                inputDS = BuildPyramidLevel(levelMosaicInfo, levelOutputTileInfo, level);
            }
        }

        private DataSource BuildPyramidLevel(MosaicInfo levelMosaicInfo, TileInfo levelOutputTileInfo, int level)
        {
            List<int> xRange = Enumerable.Range(1, levelOutputTileInfo.CountTilesX + 1).ToList();
            List<int> yRange = Enumerable.Range(1, levelOutputTileInfo.CountTilesY + 1).ToList();

            DataSource ogrds = CreateTileIndex("TileResult_" + level, _tileIndexFieldName, _sourceSRS,
                                               _tileIndexDriverTyp);

            foreach (int yIndex in yRange)
            {
                foreach (int xIndex in xRange)
                {
                    int offsetY = (yIndex - 1) * levelOutputTileInfo.TileHeight;
                    int offsetX = (xIndex - 1) * levelOutputTileInfo.TileWidth;

                    int width, height;
                    if (yIndex == levelOutputTileInfo.CountTilesY)
                    {
                        height = levelOutputTileInfo.LastTileHeight;
                    }
                    else
                    {
                        height = levelOutputTileInfo.TileHeight;
                    }

                    if (xIndex == levelOutputTileInfo.CountTilesX)
                    {
                        width = levelOutputTileInfo.LastTileWidth;
                    }
                    else
                    {
                        width = levelOutputTileInfo.TileWidth;
                    }

                    string tileName = GetTileName(levelMosaicInfo, levelOutputTileInfo, xIndex, yIndex, level);
                    CreatePyramidTile(levelMosaicInfo, offsetX, offsetY, width, height, tileName, ogrds);
                }
            }

            if (!string.IsNullOrEmpty(_tileIndexName))
            {
                // TODO: check this
                string shapeName = GetTargetDir(level) + _tileIndexName;
                CopyTileIndexToDisk(ogrds, shapeName);
            }

            if (!string.IsNullOrEmpty(_csvFileName))
            {
                string csvName = GetTargetDir(level) + _csvFileName;
                CopyTileIndexToCSV(ogrds, csvName);
            }

            return ogrds;
        }

        private void CreatePyramidTile(MosaicInfo levelMosaicInfo, int offsetX, int offsetY, int width, int height, string tileName, DataSource ogrds)
        {
            double sx = levelMosaicInfo.ScaleX * _scaleFactor;
            double sy = levelMosaicInfo.ScaleY * _scaleFactor;

            AffineTransformDecorator dec =
                new AffineTransformDecorator(new[]
                                                 {
                                                     levelMosaicInfo.Ulx + offsetX*sx, sx, 0,
                                                     levelMosaicInfo.Uly + offsetY*sy, 0, sy
                                                 });

            Dataset sFh = levelMosaicInfo.GetDataSet((int) Math.Round(dec.Ulx),
                                                     (int) Math.Round((dec.Uly + height*dec.ScaleY)),
                                                     (int) Math.Round((dec.Ulx + width*dec.ScaleX)),
                                                     (int) Math.Round(dec.Uly));

            if (sFh == null)
            {
                return;
            }

            if (ogrds != null)
            {
                var points = dec.PointsFor(width, height);
                AddFeature(ref ogrds, tileName, points);
            }

            DataType bt;
            if (_bandType == DataType.GDT_Unknown)
            {
                bt = levelMosaicInfo.BandType;
            }
            else
            {
                bt = _bandType;
            }

            double[] geotransform = { dec.Ulx, dec.ScaleX, 0, dec.Uly, 0, dec.ScaleY };
            int bands = levelMosaicInfo.Bands;

            Dataset tFh;
            if (_memDriver == null)
            {
                tFh = _driver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
            }
            else
            {
                tFh = _memDriver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
            }

            if (tFh == null)
            {
                throw new Exception("Creation failed, terminating gdal_tile.");
            }

            tFh.SetGeoTransform(geotransform);
            tFh.SetProjection(levelMosaicInfo.Projection);

            for (int band = 1; band < bands + 1; band++)
            {
                Band tBand = tFh.GetRasterBand(band);
                if (levelMosaicInfo.Ct != null)
                {
                    tBand.SetRasterColorTable(levelMosaicInfo.Ct);
                }
                tBand.SetRasterColorInterpretation(levelMosaicInfo.Ci[band - 1]);
            }

            // TODO: check this
            var res = Gdal.ReprojectImage(sFh, tFh, null, null, _resamplingMethod, 0, 0, null, null);

            if (res != 0)
            {
                throw new Exception("Reprojection failed for " + tileName + ", error " + res + ".");
            }

            levelMosaicInfo.CloseDataSet(sFh);

            if (_memDriver != null)
            {
                Dataset ttFh = _driver.CreateCopy(tileName, tFh, 0, _createOptions.ToArray(), null, null);
                ttFh.Dispose();
            }

            tFh.Dispose();
        }

        private DataSource TileImage(MosaicInfo mInfo, TileInfo ti)
        {
            _lastRowIndex = -1;
            DataSource ogrds = CreateTileIndex("TileResult_0", _tileIndexFieldName, _sourceSRS, _tileIndexDriverTyp);

            for (int yIndex = 1; yIndex < ti.CountTilesY + 1; yIndex++)
            {
                for (int xIndex = 1; xIndex < ti.CountTilesX + 1; xIndex++)
                {
                    int width, height;
                    int offsetY = (yIndex - 1) * ti.TileHeight;
                    int offsetX = (xIndex - 1) * ti.TileWidth;
                    if (yIndex == ti.CountTilesY)
                    {
                        height = ti.LastTileHeight;
                    }
                    else
                    {
                        height = ti.TileHeight;
                    }

                    if (xIndex == ti.CountTilesX)
                    {
                        width = ti.LastTileWidth;
                    }
                    else
                    {
                        width = ti.TileWidth;
                    }

                    string tileName = string.Empty;
                    if (_useDirForEachRow)
                    {
                        tileName = GetTileName(mInfo, ti, xIndex, yIndex, 0);
                    }
                    else
                    {
                        tileName = GetTileName(mInfo, ti, xIndex, yIndex);
                    }

                    CreateTile(ref mInfo, offsetX, offsetY, width, height, tileName, ref ogrds);
                }
            }

            if (!string.IsNullOrEmpty(_tileIndexName))
            {
                string shapeName = string.Empty;
                if (_useDirForEachRow && !_pyramidOnly)
                {
                    shapeName = GetTargetDir(0) + _tileIndexName;
                }
                else
                {
                    shapeName = GetTargetDir() + _tileIndexName;
                }

                CopyTileIndexToDisk(ogrds, shapeName);
            }

            if (!string.IsNullOrEmpty(_csvFileName))
            {
                string csvName = string.Empty;
                if (_useDirForEachRow && !_pyramidOnly)
                {
                    csvName = GetTargetDir(0) + _csvFileName;
                }
                else
                {
                    csvName = GetTargetDir() + _csvFileName;
                }

                CopyTileIndexToCSV(ogrds, csvName);
            }

            return ogrds;
        }

        private void CopyTileIndexToDisk(DataSource ogrds, string fileName)
        {
            DataSource shapeDS = CreateTileIndex(fileName, _tileIndexFieldName, ogrds.GetLayerByIndex(0).GetSpatialRef(),
                                                 "ESRI Shapefile");
            ogrds.GetLayerByIndex(0).ResetReading();

            while (true)
            {
                Feature feature = ogrds.GetLayerByIndex(0).GetNextFeature();
                if (feature == null)
                {
                    break;
                }

                Feature newFeature = feature.Clone();
                string basename = Path.GetFileName(feature.GetFieldAsString(0));
                if (_useDirForEachRow)
                {
                    string fExt = Path.GetExtension(basename);
                    basename = fExt + "/" + basename;
                }

                newFeature.SetField(0, basename);
                shapeDS.GetLayerByIndex(0).CreateFeature(newFeature);
            }

            CloseTileIndex(shapeDS);

            ogrds.Dispose();
        }

        private void CopyTileIndexToCSV(DataSource ogrds, string fileName)
        {
            using (File.Open(fileName, FileMode.Create))
            {
            }

            ogrds.GetLayerByIndex(0).ResetReading();

            while (true)
            {
                Feature feature = ogrds.GetLayerByIndex(0).GetNextFeature();
                if (feature == null)
                {
                    break;
                }

                string basename = Path.GetFileName(feature.GetFieldAsString(0));
                if (_useDirForEachRow)
                {
                    string fExt = Path.GetExtension(basename);
                    basename = fExt + "/" + basename;
                }

                using (StreamWriter sw = new StreamWriter(fileName, true))
                {
                    sw.Write(basename);
                    Geometry geom = feature.GetGeometryRef();
                    Envelope coords = new Envelope();
                    geom.GetEnvelope(coords);

                    // TODO: check this
                    sw.Write(_csvDelimiter + coords.MinX);
                    sw.Write(_csvDelimiter + coords.MaxX);
                    sw.Write(_csvDelimiter + coords.MinY);
                    sw.Write(_csvDelimiter + coords.MaxY);
                    sw.Write(Environment.NewLine);
                }
            }

            ogrds.Dispose();
        }

        private void CloseTileIndex(DataSource ogrDataSource)
        {
            ogrDataSource.Dispose();
        }

        private void CreateTile(ref MosaicInfo mInfo, int offsetX, int offsetY, int width, int height, string tileName, ref DataSource ogrds)
        {
            DataType bt;
            if (_bandType == DataType.GDT_Unknown)
            {
                bt = mInfo.BandType;
            }
            else
            {
                bt = _bandType;
            }

            AffineTransformDecorator dec = new AffineTransformDecorator(new[]
                                                                            {
                                                                                mInfo.Ulx,
                                                                                mInfo.ScaleX,
                                                                                0d,
                                                                                mInfo.Uly,
                                                                                0d,
                                                                                mInfo.ScaleY
                                                                            });

            Dataset sFh = mInfo.GetDataSet(
                (int) Math.Round((dec.Ulx + offsetX*dec.ScaleX)),
                (int) Math.Round((dec.Uly + offsetY*dec.ScaleY + height*dec.ScaleY)),
                (int) Math.Round((dec.Ulx + offsetX*dec.ScaleX + width*dec.ScaleX)),
                (int) Math.Round((dec.Uly + offsetY*dec.ScaleY)));

            if (sFh == null)
            {
                return;
            }

            var geotransform = new[]
                                   {
                                       dec.Ulx + offsetX * dec.ScaleX,
                                       dec.ScaleX,
                                       0d,
                                       dec.Uly + offsetY * dec.ScaleY,
                                       0d,
                                       dec.ScaleY
                                   };

            if (ogrds != null)
            {
                AffineTransformDecorator dec2 = new AffineTransformDecorator(geotransform);
                var points = dec2.PointsFor(width, height);
                AddFeature(ref ogrds, tileName, points);
            }

            int bands = mInfo.Bands;

            Dataset tFh;
            if (_memDriver == null)
            {
                tFh = _driver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
            }
            else
            {
                tFh = _memDriver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
            }

            if (tFh == null)
            {
                throw new Exception("Creation failed, terminating gdal_tile.");
            }

            tFh.SetGeoTransform(geotransform);
            if (_sourceSRS != null)
            {
                string wkt = string.Empty;
                _sourceSRS.ExportToWkt(out wkt);
                tFh.SetProjection(wkt);
            }

            int readX = Math.Min(sFh.RasterXSize, width);
            int readY = Math.Min(sFh.RasterYSize, height);

            Band tBand, sBand;
            for (int band = 1; band < bands + 1; band++)
            {
                sBand = sFh.GetRasterBand(band);
                tBand = tFh.GetRasterBand(band);

                if (mInfo.Ct != null)
                {
                    tBand.SetRasterColorTable(mInfo.Ct);
                }

                // TODO: check this
                byte[] data = new byte[readX * readY];
                sBand.ReadRaster(0, 0, readX, readY, data, readX, readY, 0, 0);
                tBand.WriteRaster(0, 0, readX, readY, data, readX, readY, 0, 0);
            }


            mInfo.CloseDataSet(sFh);

            Dataset ttFh;
            if (_memDriver != null)
            {
                ttFh = _driver.CreateCopy(tileName, tFh, 0, _createOptions.ToArray(), null, null);
                ttFh.Dispose();
            }

            tFh.Dispose();
            sFh.Dispose();
        }

        private string GetTileName(MosaicInfo mInfo, TileInfo ti, int xIndex, int yIndex, int level = -1)
        {
            int max = ti.CountTilesX;
            if (ti.CountTilesY > max)
            {
                max = ti.CountTilesY;
            }

            int countDigits = max.ToString().Length;
            string fileName = Path.GetFileNameWithoutExtension(mInfo.FileName).Replace("@", "");
            string extension = Path.GetExtension(mInfo.FileName).Replace("@", "");

            string format = string.Empty;
            if (_useDirForEachRow)
            {
                //format = Path.Combine(GetTargetDir(level) + yIndex,
                //                             fileName + "_%0" + countDigits + "i_%0" + countDigits + "i");
                format = Path.Combine(GetTargetDir(level) + yIndex,
                                             fileName + "_" + yIndex.ToString().PadLeft(countDigits, '0') + "_" + xIndex.ToString().PadLeft(countDigits, '0'));

                // see if there was a switch in the row, if so then create new dir for row
                if (_lastRowIndex < yIndex)
                {
                    _lastRowIndex = yIndex;
                    if (!Directory.Exists(GetTargetDir(level) + yIndex))
                    {
                        Directory.CreateDirectory(GetTargetDir(level) + yIndex);
                    }
                }
            }
            else
            {
                //format = GetTargetDir(level) + fileName + "_%0" + countDigits + "i_%0" + countDigits + "i";
                format = GetTargetDir(level) + fileName + "_" + yIndex.ToString().PadLeft(countDigits, '0') + "_" + xIndex.ToString().PadLeft(countDigits, '0');
            }

            // check for the extension that should be used
            if (string.IsNullOrEmpty(extension))
            {
                format += extension;
            }
            else
            {
                format += "." + _extension;
            }

            return format;
        }

        private string GetTargetDir(int level = -1)
        {
            string directory = _targetDir;
            if (!directory.EndsWith("\\"))
            {
                directory = directory + "\\";
            }
            return level == -1 ? directory : (Path.Combine(directory, level.ToString()) + "\\");
        }

        private DataSource GetTileIndexFromFiles(IEnumerable<string> inputTiles, string driverType)
        {
            DataSource ogrTileIndexDS = CreateTileIndex("TileIndex", _tileIndexFieldName, null, driverType);
            foreach (string inputTile in inputTiles)
            {
                Dataset fhInputTile = Gdal.Open(inputTile, Access.GA_ReadOnly);
                if (fhInputTile == null)
                {
                    return null;
                }

                // TODO: check this
                double[] argout = new double[6];
                fhInputTile.GetGeoTransform(argout);
                AffineTransformDecorator dec = new AffineTransformDecorator(argout);
                List<Point> points = dec.PointsFor(fhInputTile.RasterXSize, fhInputTile.RasterYSize);

                if (_rasterSize == -1 || _meters == -1)
                {
                    _meters = dec.ScaleX * fhInputTile.RasterXSize;
                    _rasterSize = fhInputTile.RasterXSize;
                }

                AddFeature(ref ogrTileIndexDS, inputTile, points);
                fhInputTile.Dispose(); // TODO: check this
            }

            if (_verbose)
            {
                Debug.WriteLine("Finished");
            }

            return ogrTileIndexDS;
        }

        private void AddFeature(ref DataSource ogrDataSource, string location, List<Point> points)
        {
            Layer ogrLayer = ogrDataSource.GetLayerByIndex(0); // TODO: check this
            Feature ogrFeature = new Feature(ogrLayer.GetLayerDefn());
            if (ogrFeature == null)
            {
                throw new Exception("Could not vreate the feature.");
            }

            ogrFeature.SetField(_tileIndexFieldName, location);

            string wkt = string.Format("POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{8} {9} ))",
                points.ElementAt(0).X,
                points.ElementAt(0).Y,
                points.ElementAt(1).X,
                points.ElementAt(1).Y,
                points.ElementAt(2).X,
                points.ElementAt(2).Y,
                points.ElementAt(3).X,
                points.ElementAt(3).Y,
                points.ElementAt(0).X,
                points.ElementAt(0).Y);

            var ogrGeometry = Ogr.CreateGeometryFromWkt(ref wkt, ogrLayer.GetSpatialRef());
            if (ogrGeometry == null)
            {
                throw new Exception("Could not create geometry.");
            }

            ogrFeature.SetGeometryDirectly(ogrGeometry);
            ogrLayer.CreateFeature(ogrFeature);
            ogrFeature.Dispose(); // TODO: check this
        }

        private static DataSource CreateTileIndex(string dsName, string fieldName, SpatialReference srs, string driverName)
        {
            var ogrDriver = Ogr.GetDriverByName(driverName);
            if (ogrDriver == null)
            {
                throw new Exception("ESRI Shapefile _driver not found.");
            }

            DataSource ogrDataSource = ogrDriver.Open(dsName, 1);
            if (ogrDataSource != null)
            {
                // python:
                // OGRDataSource.Destroy()
                ogrDataSource.Dispose(); // TODO: check this
                ogrDriver.DeleteDataSource(dsName);
            }

            ogrDataSource = ogrDriver.CreateDataSource(dsName, new string[] { }); // TODO: check this
            if (ogrDataSource == null)
            {
                throw new Exception("Could not open datasource " + dsName);
            }

            Layer ogrLayer = ogrDataSource.CreateLayer("index", srs, wkbGeometryType.wkbPolygon, new string[] { }); // TODO: check this
            if (ogrLayer == null)
            {
                throw new Exception("Could not create layer.");
            }

            FieldDefn ogrFieldDefn = new FieldDefn(fieldName, FieldType.OFTString);
            if (ogrFieldDefn == null)
            {
                throw new Exception("Could not create FieldDefn for " + fieldName);
            }

            ogrFieldDefn.SetWidth(256);
            if (ogrLayer.CreateField(ogrFieldDefn, 1) != 0) // TODO: check this
            {
                throw new Exception("Could not create field for " + fieldName);
            }

            return ogrDataSource;
        }

        private static void Usage()
        {
            Console.WriteLine("Usage: gdal_retile ");
            Console.WriteLine("        [-v] [-co NAME=VALUE]* [-of out_format]");
            Console.WriteLine("        [-ps pixelWidth pixelHeight]");
            Console.WriteLine("        [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/");
            Console.WriteLine("               CInt16/CInt32/CFloat32/CFloat64}]");
            Console.WriteLine("        [-tileIndex _tileIndexName [-tileIndexField fieldName]]");
            Console.WriteLine("        [-csv fileName [-csvDelim delimiter]]");
            Console.WriteLine("        [-s_srs srs_def]  [-pyramidOnly] -levels numberoflevels");
            Console.WriteLine("        [-r {near/bilinear/cubic/cubicspline/lanczos}]");
            Console.WriteLine("        [-useDirForEachRow]");
            Console.WriteLine("        [-targetDir tile_directory input_files]");

        }

        private void UsageFormat()
        {
            Console.WriteLine("Valid formats:");
            int count = Gdal.GetDriverCount();
            for (int i = 0; i < count; i++)
            {
                _driver = Gdal.GetDriver(i);
                Console.WriteLine(_driver.ShortName);
            }
        }

        // a class providing some usefull methods for affine Transformations
        internal class AffineTransformDecorator
        {
            private double[] _geotransform;
            private double _scaleX, _scaleY, _ulx, _uly;

            public double ScaleX
            {
                get { return _scaleX; }
                set { _scaleX = value; }
            }

            public double ScaleY
            {
                get { return _scaleY; }
                set { _scaleY = value; }
            }

            public double Ulx
            {
                get { return _ulx; }
                set { _ulx = value; }
            }

            public double Uly
            {
                get { return _uly; }
                set { _uly = value; }
            }

            public AffineTransformDecorator(double[] transform)
            {
                _geotransform = transform;
                _scaleX = _geotransform[1];
                _scaleY = _geotransform[5];

                if (_scaleY > 0)
                    _scaleY *= -1;

                _ulx = _geotransform[0];
                _uly = _geotransform[3];
            }

            public List<Point> PointsFor(int width, int height)
            {
                List<Point> points = new List<Point>();
                int w = (int)_scaleX * width;
                int h = (int)_scaleY * height;

                points.Add(new Point((int) Math.Round(_ulx), (int) Math.Round(_uly)));
                points.Add(new Point((int) Math.Round((_ulx + w)), (int) Math.Round(_uly)));
                points.Add(new Point((int) Math.Round(_ulx + w), (int) Math.Round(_uly + h)));
                points.Add(new Point((int) Math.Round(_ulx), (int) Math.Round(_uly + h)));

                return points;
            }
        }

        // a class holding information about a GDAL file or a GDAL fileset
        internal class MosaicInfo : IDisposable
        {
            private Driver _tempDriver;
            private string _fileName;
            private DataSetCache _cache;
            private DataSource _ogrTileIndexDS;
            private int _bands;
            private DataType _bandType;
            private string _projection;
            private double _scaleX, _scaleY;
            private ColorTable _ct;
            private ColorInterp[] _ci;
            private double _ulx, _uly, _lrx, _lry;
            private double _xSize, _ySize;

            public double XSize
            {
                get { return _xSize; }
                set { _xSize = value; }
            }

            public double YSize
            {
                get { return _ySize; }
                set { _ySize = value; }
            }

            public double ScaleX
            {
                get { return _scaleX; }
                set { _scaleX = value; }
            }

            public double ScaleY
            {
                get { return _scaleY; }
                set { _scaleY = value; }
            }

            public string Projection
            {
                get { return _projection; }
                set { _projection = value; }
            }

            public string FileName
            {
                get { return _fileName; }
                set { _fileName = value; }
            }

            public DataType BandType
            {
                get { return _bandType; }
                set { _bandType = value; }
            }

            public double Ulx
            {
                get { return _ulx; }
                set { _ulx = value; }
            }

            public double Uly
            {
                get { return _uly; }
                set { _uly = value; }
            }

            public double Lrx
            {
                get { return _lrx; }
                set { _lrx = value; }
            }

            public double Lry
            {
                get { return _lry; }
                set { _lry = value; }
            }

            public int Bands
            {
                get { return _bands; }
                set { _bands = value; }
            }

            public ColorTable Ct
            {
                get { return _ct; }
                set { _ct = value; }
            }

            public ColorInterp[] Ci
            {
                get { return _ci; }
                set { _ci = value; }
            }

            public MosaicInfo(string fileName, DataSource inputDS)
            {
                _tempDriver = Gdal.GetDriverByName("MEM");
                _fileName = fileName;
                _cache = new DataSetCache();
                _ogrTileIndexDS = inputDS;

                // TODO: check this
                // python: self.ogrTileIndexDS.GetLayer().ResetReading()
                _ogrTileIndexDS.GetLayerByIndex(0).ResetReading();
                Feature feature = _ogrTileIndexDS.GetLayerByIndex(0).GetNextFeature();
                string imgLocation = feature.GetFieldAsString(0);

                Dataset fhInputTile = _cache.Get(imgLocation);

                _bands = fhInputTile.RasterCount;
                _bandType = fhInputTile.GetRasterBand(1).DataType;
                _projection = fhInputTile.GetProjection();

                double[] argout = new double[6];
                fhInputTile.GetGeoTransform(argout);
                AffineTransformDecorator dec = new AffineTransformDecorator(argout);

                _scaleX = dec.ScaleX;
                _scaleY = dec.ScaleY;

                ColorTable ct = fhInputTile.GetRasterBand(1).GetRasterColorTable();
                if (ct != null)
                {
                    _ct = ct.Clone();
                }
                else
                {
                    _ct = null;
                }

                List<ColorInterp> ciList = new List<ColorInterp>();
                for (int i = 0; i < _bands; i++)
                {
                    ciList.Add(ColorInterp.GCI_Undefined);
                }

                if (_ci == null)
                {
                    _ci = new ColorInterp[_bands];
                }

                for (int iBand = 0; iBand < _bands; iBand++)
                {
                    _ci[iBand] = fhInputTile.GetRasterBand(iBand + 1).GetRasterColorInterpretation();
                }

                // TODO: check this
                // python:
                //                extent = self.ogrTileIndexDS.GetLayer().GetExtent()
                //                self.ulx = extent[0]
                //                self.uly = extent[3]
                //                self.lrx = extent[1]
                //                self.lry = extent[2]
                Envelope extent = new Envelope();
                _ogrTileIndexDS.GetLayerByIndex(0).GetExtent(extent, 1);

                _ulx = extent.MinX;
                _uly = extent.MaxY;
                _lrx = extent.MaxX;
                _lry = extent.MinY;

                _xSize = (int) Math.Round((Math.Round((_lrx - _ulx)/_scaleX)));
                _ySize = Math.Abs((int) Math.Round((Math.Round((_uly - _lry)/_scaleY))));
            }

            public Dataset GetDataSet(double minX, double minY, double maxX, double maxY)
            {
                _ogrTileIndexDS.GetLayerByIndex(0).ResetReading();
                _ogrTileIndexDS.GetLayerByIndex(0).SetSpatialFilterRect(minX, minY, maxX, maxY);

                List<Feature> features = new List<Feature>();
                Envelope envelope = null;

                while (true)
                {
                    Feature feature = _ogrTileIndexDS.GetLayerByIndex(0).GetNextFeature();
                    if (feature == null)
                    {
                        break;
                    }
                    features.Add(feature);
                    if (envelope == null)
                    {
                        // TODO: check this
                        // python:
                        // envelope=feature.GetGeometryRef().GetEnvelope()
                        envelope = new Envelope();
                        feature.GetGeometryRef().GetEnvelope(envelope);
                    }
                    else
                    {
                        Envelope featureEnv = new Envelope();
                        feature.GetGeometryRef().GetEnvelope(featureEnv);

                        envelope.MinX = Math.Min(featureEnv.MinX, envelope.MinX);
                        envelope.MinY = Math.Min(featureEnv.MinY, envelope.MinY);
                        envelope.MaxX = Math.Max(featureEnv.MaxX, envelope.MaxX);
                        envelope.MaxY = Math.Max(featureEnv.MaxY, envelope.MaxY);
                    }
                }

                if (envelope == null)
                {
                    return null;
                }

                // enlarge to query rect if necessairy
                envelope.MinX = Math.Min(minX, envelope.MinX);
                envelope.MinY = Math.Min(minY, envelope.MinY);
                envelope.MaxX = Math.Max(maxX, envelope.MaxX);
                envelope.MaxY = Math.Max(maxY, envelope.MaxY);

                _ogrTileIndexDS.GetLayerByIndex(0).SetSpatialFilter(null);

                // merge tiles
                int resultSizeX = (int) Math.Round((Math.Ceiling((double) (maxX - minX)/(double) _scaleX)));
                int resultSizeY = (int) Math.Round((Math.Ceiling((double) (minY - maxY)/(double) _scaleY)));

                Dataset resultDS = _tempDriver.Create("TEMP", resultSizeX, resultSizeY, _bands, _bandType, new string[] { });
                resultDS.SetGeoTransform(new[] { minX, _scaleX, 0, maxY, 0, _scaleY });

                foreach (Feature feature in features)
                {
                    string featureName = feature.GetFieldAsString(0);
                    Dataset sourceDS = _cache.Get(featureName);

                    double[] argout = new double[6];
                    sourceDS.GetGeoTransform(argout);
                    AffineTransformDecorator dec = new AffineTransformDecorator(argout);

                    // calculate read and write offsets
                    int readOffsetX = (int) Math.Round((Math.Round((double) (minX - dec.Ulx)/(double) _scaleX)));
                    int readOffsetY = (int) Math.Round((Math.Round((double) (maxY - dec.Uly)/(double) _scaleY)));

                    int writeOffsetX = 0;
                    if (readOffsetX < 0)
                    {
                        writeOffsetX = readOffsetX * (-1);
                        readOffsetX = 0;
                    }
                    int writeOffsetY = 0;
                    if (readOffsetY < 0)
                    {
                        writeOffsetY = readOffsetY * (-1);
                        readOffsetY = 0;
                    }

                    // calculate read and write dimensions
                    int readX = Math.Min(Math.Min(resultSizeX, sourceDS.RasterXSize - readOffsetX), resultSizeX - writeOffsetX);
                    if (readX <= 0)
                        continue;

                    int readY = Math.Min(Math.Min(resultSizeY, sourceDS.RasterYSize - readOffsetY), resultSizeY - writeOffsetY);
                    if (readY <= 0)
                        continue;

                    for (int bandNr = 1; bandNr < _bands + 1; bandNr++)
                    {
                        Band sBand = sourceDS.GetRasterBand(bandNr);
                        Band tBand = resultDS.GetRasterBand(bandNr);

                        if (_ct != null)
                        {
                            tBand.SetRasterColorTable(_ct);
                        }
                        tBand.SetRasterColorInterpretation(_ci[bandNr - 1]);
                        // TODO: check this
                        //var data = sBand.ReadRaster(readOffsetX, readOffsetY, readX, readY, readX, readY, _bandType);
                        byte[] buffer = new byte[readX * readY];
                        sBand.ReadRaster(readOffsetX, readOffsetY, readX, readY, buffer, readX, readY, 0, 0);
                        tBand.WriteRaster(writeOffsetX, writeOffsetY, readX, readY, buffer, readX, readY, 0, 0);
                    }
                }

                return resultDS;
            }

            public void CloseDataSet(Dataset memDS)
            {
                memDS.Dispose();
            }

            public void Dispose()
            {
                _cache.Dispose();
                _ogrTileIndexDS.Dispose();
            }
        }

        internal class DataSetCache : IDisposable
        {
            private int _cacheSize;
            private Queue<string> _queue;
            private Dictionary<string, Dataset> _dict;

            public DataSetCache()
            {
                _cacheSize = 8;
                _queue = new Queue<string>();
                _dict = new Dictionary<string, Dataset>();
            }

            public Dataset Get(string name)
            {
                if (_dict.ContainsKey(name))
                    return _dict[name];

                Dataset result = Gdal.Open(name, Access.GA_ReadOnly);
                if (result == null)
                {
                    throw new Exception("Error opening: " + name);
                }

                if (_queue.Count() == _cacheSize)
                {
                    string toRemove = _queue.Dequeue();
                    _dict.Remove(toRemove);
                }

                _queue.Enqueue(name);
                _dict[name] = result;

                return result;
            }

            public void Dispose()
            {
                _dict.Clear();
                _queue.Clear();
                _dict = null;
                _queue = null;
                GC.Collect();
            }
        }

        // a class holding info how to tile
        internal class TileInfo
        {
            private int _tileWidth, _tileHeight;
            private int _countTilesX, _countTilesY;
            private int _lastTileWidth, _lastTileHeight;

            public int CountTilesX
            {
                get { return _countTilesX; }
                set { _countTilesX = value; }
            }

            public int CountTilesY
            {
                get { return _countTilesY; }
                set { _countTilesY = value; }
            }

            public int TileWidth
            {
                get { return _tileWidth; }
                set { _tileWidth = value; }
            }

            public int TileHeight
            {
                get { return _tileHeight; }
                set { _tileHeight = value; }
            }

            public int LastTileWidth
            {
                get { return _lastTileWidth; }
                set { _lastTileWidth = value; }
            }

            public int LastTileHeight
            {
                get { return _lastTileHeight; }
                set { _lastTileHeight = value; }
            }

            public TileInfo(double xSize, double ySize, int tileWidth, int tileHeight)
            {
                _tileWidth = tileWidth;
                _tileHeight = tileHeight;
                _countTilesX = (int) Math.Round(xSize/(double) tileWidth);
                _countTilesY = (int) Math.Round(ySize/(double) tileHeight);
                _lastTileWidth = (int) Math.Round(xSize - (double) _countTilesX*(double) tileWidth);
                _lastTileHeight = (int) Math.Round(ySize - (double) _countTilesY*(double) tileHeight);

                if (_lastTileWidth > 0)
                {
                    _countTilesX = _countTilesX + 1;
                }
                else
                {
                    _lastTileWidth = tileWidth;
                }

                if (_lastTileHeight > 0)
                {
                    _countTilesY = _countTilesY + 1;
                }
                else
                {
                    _lastTileHeight = tileHeight;
                }
            }

            public void Report()
            {
                Console.WriteLine("tileWidth       " + _tileWidth);
                Console.WriteLine("tileHeight      " + _tileHeight);
                Console.WriteLine("countTilesX:    " + _countTilesX);
                Console.WriteLine("countTilesY:    " + _countTilesY);
                Console.WriteLine("lastTileWidth:  " + _lastTileWidth);
                Console.WriteLine("lastTileHeight: " + _lastTileHeight);
            }
        }
    }
}