Advertisement
Guest User

Untitled

a guest
Oct 8th, 2012
223
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C# 76.81 KB | None | 0 0
  1. using System;
  2. using System.Collections.Generic;
  3. using System.Diagnostics;
  4. using System.Drawing;
  5. using System.IO;
  6. using System.Linq;
  7. using System.Text.RegularExpressions;
  8. using System.Windows;
  9. using System.Windows.Media.Imaging;
  10. using System.Text;
  11. using System.Xml.Linq;
  12. using OSGeo.GDAL;
  13. using OSGeo.OGR;
  14. using OSGeo.OSR;
  15. using Driver = OSGeo.GDAL.Driver;
  16. using Point = System.Drawing.Point;
  17. using Size = System.Drawing.Size;
  18.  
  19. namespace RasterLoader
  20. {
  21.     public class GdalRetile
  22.     {
  23.         private const double MetersPerInch = 0.0254;
  24.         private const double DPI = 96;
  25.  
  26.         private string _format = string.Empty;
  27.         private List<string> _createOptions = new List<string>();
  28.         private DataType _bandType;
  29.         private bool _verbose;
  30.         private string _targetDir = string.Empty;
  31.         private int _tileWidth = -1;
  32.         private int _tileHeight = -1;
  33.         private string _resamplingMethodString = string.Empty;
  34.         private ResampleAlg _resamplingMethod;
  35.         private int _levels = -1;
  36.         private SpatialReference _sourceSRS;
  37.         private string _tileIndexName = string.Empty;
  38.         private string _tileIndexFieldName = string.Empty;
  39.         private string _tileIndexDriverTyp = string.Empty;
  40.         private string _csvFileName = string.Empty;
  41.         private string _csvDelimiter = string.Empty;
  42.         private List<string> _names = new List<string>();
  43.         private Driver _driver;
  44.         private Driver _memDriver;
  45.         private int _lastRowIndex;
  46.         private string _extension;
  47.         private bool _8To32Bit;
  48.         private double _scaleFactor;
  49.         private double _rasterSize = -1;
  50.         private double _meters = -1;
  51.         private string _layerFileName;
  52.         private bool _pyramidOnly;
  53.         private bool _useDirForEachRow;
  54.  
  55.         /// <summary>
  56.         /// Creation options for output file. Multiple options can be specified.
  57.         /// </summary>
  58.         public enum CreationOptions
  59.         {
  60.             /// <summary>
  61.             /// By default stripped TIFF files are created. This option can be used to force creation of tiled TIFF files.
  62.             /// </summary>
  63.             TILED,
  64.             /// <summary>
  65.             /// TIFF files can contain internal transparency masks. The GeoTIFF driver recognizes an internal directory as being a transparency mask when the FILETYPE_MASK bit value is set on the TIFFTAG_SUBFILETYPE tag. According to the TIFF specification, such internal transparency masks contain 1 sample of 1-bit data. Although the TIFF specification allows for higher resolutions for the internal transparency mask, the GeoTIFF driver only supports internal transparency masks of the same dimensions as the main image. Transparency masks of internal overviews are also supported.
  66.             /// </summary>
  67.             GDAL_TIFF_INTERNAL_MASK,
  68.             /// <summary>
  69.             /// Set the compression to use: JPEG (JPEG files).
  70.             /// </summary>
  71.             COMPRESS_JPEG,
  72.             /// <summary>
  73.             /// Set the compression to use: LZW (PNG files).
  74.             /// </summary>
  75.             COMPRESS_LZW,
  76.             /// <summary>
  77.             /// Set the JPEG quality [1 - 100].
  78.             /// </summary>
  79.             JPEG_QUALITY,
  80.             /// <summary>
  81.             /// Set the JPEG quality [1 - 100].
  82.             /// </summary>
  83.             JPEG_QUALITY_OVERVIEW,
  84.             /// <summary>
  85.             /// The first "extrasample" is marked as being alpha if there are any extra samples. This is necessary if you want to produce a greyscale TIFF file with an alpha band (for instance).
  86.             /// </summary>
  87.             ALPHA
  88.         }
  89.  
  90.         /// <summary>
  91.         /// Output format, defaults to GeoTIFF (GTiff).
  92.         /// </summary>
  93.         public string OutputFormat
  94.         {
  95.             get { return _format; }
  96.             set { _format = value; }
  97.         }
  98.  
  99.         /// <summary>
  100.         /// Force the output image bands to have a specific type. Use type names (ie. Byte, Int16,...).
  101.         /// </summary>
  102.         public string BandDataType
  103.         {
  104.             get { return Gdal.GetDataTypeName(_bandType); }
  105.             set
  106.             {
  107.                 _bandType = Gdal.GetDataTypeByName(value);
  108.                 if (_bandType == DataType.GDT_Unknown)
  109.                 {
  110.                     throw new Exception("Unknown GDAL data type: " + value);
  111.                 }
  112.             }
  113.         }
  114.  
  115.         /// <summary>
  116.         /// The directory where the tile result is created. Pyramids are stored in subdirs numbered from 1. Created tile names have a numbering schema and contain the name of the source tiles(s).
  117.         /// </summary>
  118.         public string TargetDirectory
  119.         {
  120.             get { return _targetDir; }
  121.             set
  122.             {
  123.                 _targetDir = value;
  124.                 if (!Directory.Exists(_targetDir))
  125.                 {
  126.                     throw new Exception("TargetDir " + _targetDir + " does not exist.");
  127.                 }
  128.             }
  129.         }
  130.  
  131.         /// <summary>
  132.         /// Pixel size to be used for the output file. If not specified, 512 x 512 is the default.
  133.         /// </summary>
  134.         public Size TileSize
  135.         {
  136.             get { return new Size(_tileWidth, _tileHeight); }
  137.             set
  138.             {
  139.                 _tileWidth = value.Width;
  140.                 _tileHeight = value.Height;
  141.             }
  142.         }
  143.  
  144.         /// <summary>
  145.         /// Resampling algorithm, default is nearest-neighbour.
  146.         /// </summary>
  147.         public string ResamplingAlgorithm
  148.         {
  149.             get
  150.             {
  151.                 switch (_resamplingMethod)
  152.                 {
  153.                     case ResampleAlg.GRA_NearestNeighbour:
  154.                         return "near";
  155.                     case ResampleAlg.GRA_Bilinear:
  156.                         return "bilinear";
  157.                     case ResampleAlg.GRA_Cubic:
  158.                         return "cubic";
  159.                     case ResampleAlg.GRA_CubicSpline:
  160.                         return "cubicspline";
  161.                     case (ResampleAlg)4:
  162.                         return "lanczos";
  163.                     default:
  164.                         throw new ArgumentOutOfRangeException();
  165.                 }
  166.             }
  167.             set
  168.             {
  169.                 switch (value)
  170.                 {
  171.                     case "near":
  172.                         _resamplingMethod = ResampleAlg.GRA_NearestNeighbour;
  173.                         break;
  174.                     case "bilinear":
  175.                         _resamplingMethod = ResampleAlg.GRA_Bilinear;
  176.                         break;
  177.                     case "cubic":
  178.                         _resamplingMethod = ResampleAlg.GRA_Cubic;
  179.                         break;
  180.                     case "cubicspline":
  181.                         _resamplingMethod = ResampleAlg.GRA_CubicSpline;
  182.                         break;
  183.                     case "lanczos":
  184.                         _resamplingMethod = (ResampleAlg) 4;
  185.                         break;
  186.                     default:
  187.                         throw new ArgumentOutOfRangeException(value);
  188.                 }
  189.             }
  190.         }
  191.  
  192.         /// <summary>
  193.         /// Number of pyramids levels to build (if not set, algorithm finishes, when a level consist of 1 tile only).
  194.         /// </summary>
  195.         public int Levels
  196.         {
  197.             get { return _levels; }
  198.             set { _levels = value; }
  199.         }
  200.  
  201.         /// <summary>
  202.         /// Source spatial reference to use. The coordinate systems that can be passed are anything supported by the OGRSpatialReference.SetFro‐mUserInput() call, which includes EPSG PCS and GCSes (ie.EPSG:4296), PROJ.4 declarations (as above), or the name of a .prf file containing well known text. If no srs_def is given, the srs_def of the source tiles is used (if there is any). The srs_def will be propageted to created tiles (if possible) and to the optional shape file(s).
  203.         /// </summary>
  204.         public string SpatialReference
  205.         {
  206.             get
  207.             {
  208.                 string wkt = string.Empty;
  209.                 _sourceSRS.ExportToWkt(out wkt);
  210.                 return wkt;
  211.             }
  212.             set
  213.             {
  214.                 string wkt = string.Empty;
  215.                 if (Osr.GetUserInputAsWKT(value, out wkt) != 0)
  216.                 {
  217.                     throw new Exception("Invalid -s_srs: " + value);
  218.                 }
  219.                 _sourceSRS = new SpatialReference(wkt);
  220.             }
  221.         }
  222.  
  223.         /// <summary>
  224.         /// No retiling, build only the pyramids.
  225.         /// </summary>
  226.         public bool PyramidOnly
  227.         {
  228.             get { return _pyramidOnly; }
  229.             set { _pyramidOnly = value; }
  230.         }
  231.  
  232.         /// <summary>
  233.         /// Normally the tiles of the base image are stored as described in -targetDir. For large images, some file systems have performance problems if the number of files in a directory is to big, causing gdal_retile not to finish in reasonable time. Using this parameter creates a different output structure. The tiles of the base image are stored in a subdirectory called 0, the pyramids in subdirectories numbered 1,2,.... Within each of these directories another level of subdirectories is created, numbered from 0...n, depending of how many tile rows are needed for each level. Finally, a directory contains only the the tiles for one row for a specific level. For large images a performance improvement of a factor N could be achieved.
  234.         /// </summary>
  235.         public bool UseDirForEachRow
  236.         {
  237.             get { return _useDirForEachRow; }
  238.             set { _useDirForEachRow = value; }
  239.         }
  240.  
  241.         /// <summary>
  242.         /// Convert 8 bit source image into 24 bit (recommended).
  243.         /// </summary>
  244.         public bool Convert8To24Bit
  245.         {
  246.             get { return _8To32Bit; }
  247.             set { _8To32Bit = value; }
  248.         }
  249.  
  250.         public string LayerFileName
  251.         {
  252.             get { return _layerFileName; }
  253.             set { _layerFileName = value; }
  254.         }
  255.  
  256.         /// <summary>
  257.         /// Sets or removes chosen creation option.
  258.         /// </summary>
  259.         /// <param name="co">Chosen creation option.</param>
  260.         /// <param name="add">If true, adds a new creation option. If false, removes it.</param>
  261.         /// <param name="additionalValue">If a value needs to be set manually, you need to specify it (e.g. JPEG_QUALITY=50, where 50 is your additional value).</param>
  262.         public void SetCreationOption(CreationOptions co, bool add, string additionalValue = "")
  263.         {
  264.             string option = string.Empty;
  265.             switch (co)
  266.             {
  267.                 case CreationOptions.TILED:
  268.                     option = "TILED=YES";
  269.                     break;
  270.                 case CreationOptions.GDAL_TIFF_INTERNAL_MASK:
  271.                     option = "GDAL_TIFF_INTERNAL_MASK=TRUE";
  272.                     break;
  273.                 case CreationOptions.COMPRESS_JPEG:
  274.                     option = "COMPRESS=JPEG";
  275.                     break;
  276.                 case CreationOptions.COMPRESS_LZW:
  277.                     option = "COMPRESS=LZW";
  278.                     break;
  279.                 case CreationOptions.JPEG_QUALITY:
  280.                     option = "JPEG_QUALITY=" + additionalValue;
  281.                     break;
  282.                 case CreationOptions.JPEG_QUALITY_OVERVIEW:
  283.                     option = "JPEG_QUALITY_OVERVIEW=" + additionalValue;
  284.                     break;
  285.                 case CreationOptions.ALPHA:
  286.                     option = "ALPHA=YES";
  287.                     break;
  288.                 default:
  289.                     throw new ArgumentOutOfRangeException("Unknown creation option.");
  290.             }
  291.  
  292.             if (add)
  293.                 _createOptions.Add(option);
  294.             else
  295.                 _createOptions.Remove(option);
  296.         }
  297.  
  298.         /// <summary>
  299.         /// Sets input files for GdalRetile class.
  300.         /// </summary>
  301.         /// <param name="source">Source files for GdalRetile (string variable with paths - each one in a new row, or a path to the file with paths).</param>
  302.         /// <param name="isSourcePath">Specifies if the source is a path.</param>
  303.         public void SetInputFiles(string source, bool isSourcePath = false)
  304.         {
  305.             string tmp = isSourcePath ? (new StreamReader(source)).ReadToEnd() : source;
  306.             _names = tmp.Split(new[] { "\r\n", "\n" }, StringSplitOptions.RemoveEmptyEntries).ToList();
  307.         }
  308.  
  309.         public GdalRetile()
  310.         {
  311.             Gdal.AllRegister();
  312.             InitFields();
  313.         }
  314.  
  315.         public GdalRetile(string[] args)
  316.         {
  317.             Gdal.AllRegister();
  318.             InitFields();
  319.            
  320.             if (args.Any())
  321.             {
  322.                 //string[] argv = Gdal.GeneralCmdLineProcessor(args, 0);
  323.                 string[] argv = args.Clone() as string[];
  324.                 if (!argv.Any())
  325.                     throw new Exception("Command line has no parameters.");
  326.  
  327.                 int i = 1;
  328.                 while (i < argv.Count())
  329.                 {
  330.                     string arg = argv.ElementAt(i);
  331.                     string ext;
  332.                     switch (arg)
  333.                     {
  334.                         case "-of":
  335.                             i++;
  336.                             _format = argv.ElementAt(i);
  337.                             break;
  338.                         case "-ot":
  339.                             i++;
  340.                             _bandType = Gdal.GetDataTypeByName(argv.ElementAt(i));
  341.                             if (_bandType == DataType.GDT_Unknown)
  342.                             {
  343.                                 throw new Exception("Unknown GDAL data type: " + argv.ElementAt(i));
  344.                             }
  345.                             break;
  346.                         case "-co":
  347.                             i++;
  348.                             _createOptions.Add(argv.ElementAt(i));
  349.                             break;
  350.                         case "-v":
  351.                             _verbose = true;
  352.                             break;
  353.                         case "-targetDir":
  354.                             i++;
  355.                             _targetDir = argv.ElementAt(i);
  356.                             if (!Directory.Exists(_targetDir))
  357.                             {
  358.                                 throw new Exception("TargetDir " + _targetDir + " does not exist.");
  359.                             }
  360.                             break;
  361.                         case "-ps":
  362.                             i++;
  363.                             _tileWidth = int.Parse(argv.ElementAt(i));
  364.                             i++;
  365.                             _tileHeight = int.Parse(argv.ElementAt(i));
  366.                             break;
  367.                         case "-r":
  368.                             i++;
  369.                             _resamplingMethodString = argv.ElementAt(i);
  370.                             switch (_resamplingMethodString)
  371.                             {
  372.                                 case "near":
  373.                                     _resamplingMethod = ResampleAlg.GRA_NearestNeighbour;
  374.                                     break;
  375.                                 case "bilinear":
  376.                                     _resamplingMethod = ResampleAlg.GRA_Bilinear;
  377.                                     break;
  378.                                 case "cubic":
  379.                                     _resamplingMethod = ResampleAlg.GRA_Cubic;
  380.                                     break;
  381.                                 case "cubicspline":
  382.                                     _resamplingMethod = ResampleAlg.GRA_CubicSpline;
  383.                                     break;
  384.                                 case "lanczos":
  385.                                     _resamplingMethod = (ResampleAlg) 4;
  386.                                     break;
  387.                                 default:
  388.                                     throw new Exception("Unknown resampling method: " + _resamplingMethodString);
  389.                             }
  390.                             break;
  391.                         case "-levels":
  392.                             i++;
  393.                             _levels = int.Parse(argv.ElementAt(i));
  394.                             if (_levels < 1)
  395.                             {
  396.                                 throw new Exception("Invalid number of _levels: " + _levels);
  397.                             }
  398.                             break;
  399.                         case "-s_srs":
  400.                             i++;
  401.                             string wkt = string.Empty;
  402.                             if (Osr.GetUserInputAsWKT(argv.ElementAt(i), out wkt) != 0)
  403.                             {
  404.                                 throw new Exception("Invalid -s_srs: " + argv.ElementAt(i));
  405.                             }
  406.                             _sourceSRS = new SpatialReference(wkt); // TODO: check this!
  407.                             break;
  408.                         case "-pyramidOnly":
  409.                             PyramidOnly = true;
  410.                             break;
  411.                         case "-tileIndex":
  412.                             i++;
  413.                             _tileIndexName = argv.ElementAt(i);
  414.                             ext = Path.GetExtension(_tileIndexName);
  415.                             if (ext.Length == 0)
  416.                             {
  417.                                 _tileIndexName += ".shp";
  418.                             }
  419.                             break;
  420.                         case "-tileIndexField":
  421.                             i++;
  422.                             _tileIndexFieldName = argv.ElementAt(i);
  423.                             break;
  424.                         case "-csv":
  425.                             i++;
  426.                             _csvFileName = argv.ElementAt(i);
  427.                             ext = Path.GetExtension(_csvFileName);
  428.                             if (ext.Length == 0)
  429.                             {
  430.                                 _csvFileName += ".csv";
  431.                             }
  432.                             break;
  433.                         case "-csvDelim":
  434.                             i++;
  435.                             _csvDelimiter = argv.ElementAt(i);
  436.                             break;
  437.                         case "-useDirForEachRow":
  438.                             UseDirForEachRow = true;
  439.                             break;
  440.                         case "-8to32bit":
  441.                             _8To32Bit = true;
  442.                             break;
  443.                         case "-layerFileName":
  444.                             i++;
  445.                             _layerFileName = argv.ElementAt(i);
  446.                             break;
  447.                         default:
  448.                             if (arg.StartsWith("-") && !arg.Equals("--optfile"))
  449.                             {
  450.                                 Usage();
  451.                                 throw new Exception("Unrecognized command option: " + arg);
  452.                             }
  453.                             if (arg.Equals("--optfile"))
  454.                             {
  455.                                 i++;
  456.                                 using (StreamReader sr = new StreamReader(argv.ElementAt(i)))
  457.                                 {
  458.                                     _names = sr.ReadToEnd().Split(new[] {"\r\n", "\n"}, StringSplitOptions.RemoveEmptyEntries).ToList();
  459.                                 }
  460.                             }
  461.                             else
  462.                             {
  463.                                 _names.Add(arg);
  464.                             }
  465.                             break;
  466.                     }
  467.                     i++;
  468.                 }
  469.             }
  470.             else
  471.             {
  472.                 throw new Exception("Command line has no parameters.");
  473.             }
  474.         }
  475.  
  476.         private void ValidateInputData()
  477.         {
  478.             if (!_names.Any())
  479.             {
  480.                 Usage();
  481.                 throw new Exception("No input files selected.");
  482.             }
  483.  
  484.             if (_tileWidth <= 0 || _tileHeight <= 0)
  485.             {
  486.                 Usage();
  487.                 throw new Exception("Invalid tile dimension: " + _tileWidth + "x" + _tileHeight + ".");
  488.             }
  489.  
  490.             if (string.IsNullOrEmpty(_targetDir) || string.IsNullOrWhiteSpace(_targetDir))
  491.             {
  492.                 Usage();
  493.                 throw new Exception("Missing directory for tiles [-targetDir].");
  494.             }
  495.  
  496.             if (string.IsNullOrEmpty(_layerFileName))
  497.             {
  498.                 Usage();
  499.                 throw new Exception("Missing layer file name [-layerFileName].");
  500.             }
  501.         }
  502.  
  503.         public void MakePyramid(Action<int> progressCallback)
  504.         {
  505.             ValidateInputData();
  506.  
  507.             if (_format != "GTiff")
  508.             {
  509.                 _memDriver = Gdal.GetDriverByName("MEM");
  510.             }
  511.  
  512.             _driver = Gdal.GetDriverByName(_format);
  513.             if (_driver == null)
  514.             {
  515.                 UsageFormat();
  516.                 throw new Exception("Format _driver " + _format + " not found, pick a supported _driver.");
  517.             }
  518.  
  519.             string[] driverMD = _driver.GetMetadata(null); // null for the default domain
  520.             _extension = _driver.GetMetadataItem("DMD_EXTENSION", null);
  521.  
  522.             if (driverMD.Contains("DCAP_CREATE"))
  523.             {
  524.                 _memDriver = Gdal.GetDriverByName("MEM");
  525.             }
  526.  
  527.             DataSource tileIndexDS = GetTileIndexFromFiles(_names, _tileIndexDriverTyp, false);
  528.             if (tileIndexDS == null)
  529.             {
  530.                 throw new Exception("Error building tile index.");
  531.             }
  532.  
  533.             MosaicInfo mInfo = new MosaicInfo(_names[0], tileIndexDS, false);
  534.             TileInfo ti = new TileInfo(mInfo.XSize, mInfo.YSize, _tileWidth, _tileHeight);
  535.  
  536.             if (_sourceSRS == null && mInfo.Projection.Length > 0)
  537.             {
  538.                 // TODO: check this
  539.                 string wkt = string.Empty;
  540.                 Osr.GetUserInputAsWKT(mInfo.Projection, out wkt);
  541.                 _sourceSRS = new SpatialReference(wkt);
  542.  
  543.                 if (_sourceSRS.SetFromUserInput(mInfo.Projection) != 0)
  544.                 {
  545.                     throw new Exception("Invalid projection: " + mInfo.Projection);
  546.                 }
  547.             }
  548.  
  549.             if (_verbose)
  550.             {
  551.                 ti.Report();
  552.             }
  553.  
  554.             DataSource dsCreatedTileIndex;
  555.             if (!PyramidOnly)
  556.             {
  557.                 dsCreatedTileIndex = TileImage(mInfo, ti);
  558.                 tileIndexDS.Dispose();
  559.             }
  560.             else
  561.             {
  562.                 dsCreatedTileIndex = tileIndexDS;
  563.             }
  564.  
  565.             if (_levels > 0)
  566.             {
  567.                 Action<int> callback;
  568.                 if (progressCallback != null)
  569.                 {
  570.                     callback = progressCallback;
  571.                 }
  572.                 else
  573.                 {
  574.                     callback = null;
  575.                 }
  576.                 BuildPyramid(dsCreatedTileIndex, _tileWidth, _tileHeight, callback);
  577.             }
  578.         }
  579.  
  580.         private void InitFields()
  581.         {
  582.             _verbose = false;
  583.             _createOptions = new List<string>();
  584.             _names = new List<string>();
  585.             _tileWidth = 512;
  586.             _tileHeight = 512;
  587.             _format = "GTiff";
  588.             _bandType = DataType.GDT_Unknown;
  589.             _driver = null;
  590.             _extension = null;
  591.             _memDriver = null;
  592.             _tileIndexFieldName = "location";
  593.             _tileIndexName = null;
  594.             _tileIndexDriverTyp = "Memory";
  595.             _csvDelimiter = ";";
  596.             _csvFileName = null;
  597.             _sourceSRS = null;
  598.             _targetDir = null;
  599.             _resamplingMethod = ResampleAlg.GRA_NearestNeighbour;
  600.             _levels = 10000;
  601.             PyramidOnly = false;
  602.             _lastRowIndex = -1;
  603.             UseDirForEachRow = false;
  604.         }
  605.  
  606.         private static double GetScale(int meters, int pixels)
  607.         {
  608.             return meters / pixels / MetersPerInch * DPI;
  609.         }
  610.  
  611.         private void BuildPyramid(DataSource createdTileIndexDS, int tileWidth, int tileHeight, Action<int> progressCallback)
  612.         {
  613.             DataSource inputDS = createdTileIndexDS;
  614.             for (int level = 0; level < _levels + 1; level++)
  615.             {
  616.                 _lastRowIndex = -1;
  617.  
  618.                 MosaicInfo levelMosaicInfo;
  619.                 switch (level)
  620.                 {
  621.                     case 0:
  622.                         _scaleFactor = 1;
  623.                         bool isJPEGCompression = (_createOptions.Any(option => option.Contains("COMPRESS=JPEG")))
  624.                                                      ? true
  625.                                                      : false;
  626.                         levelMosaicInfo = new MosaicInfo(_layerFileName, inputDS, isJPEGCompression); // original data for level 0 (for PNG)
  627.                         break;
  628.                     case 1:
  629.                         double scale = GetScale((int) Math.Round(_meters), (int) Math.Round(_rasterSize));
  630.                         double power = Math.Log(scale, 2);
  631.                         if (power - (double) ((int) power) > 0d)
  632.                         {
  633.                             power = Math.Ceiling(power);
  634.                         }
  635.                         _scaleFactor = Math.Pow(2, power) / scale;
  636.                         levelMosaicInfo = new MosaicInfo(_layerFileName, inputDS, _8To32Bit);
  637.                         break;
  638.                     default:
  639.                         _scaleFactor = 2;
  640.                         levelMosaicInfo = new MosaicInfo(_layerFileName, inputDS, _8To32Bit);
  641.                         break;
  642.                 }
  643.  
  644.                 Action<int> callback;
  645.                 if (progressCallback != null)
  646.                 {
  647.                     callback = progressCallback;
  648.                 }
  649.                 else
  650.                 {
  651.                     callback = null;
  652.                 }
  653.                 TileInfo levelOutputTileInfo = new TileInfo(levelMosaicInfo.XSize/_scaleFactor,
  654.                                                             levelMosaicInfo.YSize/_scaleFactor, tileWidth, tileHeight);
  655.                 int minTileSizeX = 0, minTileSizeY = 0;
  656.                 string extension = string.Empty;
  657.                 inputDS = BuildPyramidLevel(levelMosaicInfo, levelOutputTileInfo, level, callback, out minTileSizeX, out minTileSizeY, out extension);
  658.                
  659.                 StoreLevelData(level, inputDS, levelOutputTileInfo.CountTilesX, levelOutputTileInfo.CountTilesY, minTileSizeX, minTileSizeY, extension);
  660.  
  661.                 levelMosaicInfo.Dispose();
  662.  
  663.                 if (levelOutputTileInfo.CountTilesX * levelOutputTileInfo.CountTilesY == 1)
  664.                 {
  665.                     // remove empty dirs
  666.                     RemoveEmptyDirectories(Path.Combine(_targetDir, Path.GetFileNameWithoutExtension(_layerFileName)));
  667.                     return;
  668.                 }
  669.             }
  670.         }
  671.  
  672.         private void StoreLevelData(int nr, DataSource inputDS, int tilesCountX, int tilesCountY, int minTileSizeX, int minTileSizeY, string extension)
  673.         {
  674.             string xmlPath = Path.Combine(_targetDir, Path.GetFileName(_layerFileName));
  675.             XDocument doc;
  676.  
  677.             if (nr == 0)
  678.             {
  679.                 doc = new XDocument();
  680.                 XElement root = new XElement("SMPyramid");
  681.                 root.SetAttributeValue("LayerName", Path.GetFileNameWithoutExtension(_layerFileName));
  682.                 root.SetAttributeValue("LayerPath", Path.Combine(_targetDir, Path.GetFileNameWithoutExtension(_layerFileName)));
  683.                 root.SetAttributeValue("Extension", extension);
  684.                 doc.Add(root);
  685.             }
  686.             else
  687.             {
  688.                 doc = XDocument.Load(xmlPath);
  689.             }
  690.            
  691.             inputDS.GetLayerByIndex(0).ResetReading();
  692.             Feature feature = inputDS.GetLayerByIndex(0).GetNextFeature();
  693.             string imgLocation = feature.GetFieldAsString(0);
  694.  
  695.             Dataset fhInputTile = Gdal.Open(imgLocation, Access.GA_ReadOnly);
  696.  
  697.             double ulx = 0, uly = 0, lrx = 0, lry = 0;
  698.             GdalHelper.GetGeocodedInfo(imgLocation, out ulx, out lrx, out lry, out uly);
  699.  
  700.             double[] argout = new double[6];
  701.             fhInputTile.GetGeoTransform(argout);
  702.  
  703.             AffineTransformDecorator dec = new AffineTransformDecorator(argout);
  704.          
  705. //            Dataset lastTile = Gdal.Open(lastTileNameForLevel, Access.GA_ReadOnly);
  706.  
  707.             XElement newNode = new XElement("Level");
  708.             newNode.SetAttributeValue("NR", nr);
  709.             newNode.SetAttributeValue("RasterXSize", fhInputTile.RasterXSize);
  710.             newNode.SetAttributeValue("RasterYSize", fhInputTile.RasterYSize);
  711.             newNode.SetAttributeValue("LastTileXSize", minTileSizeX);
  712.             newNode.SetAttributeValue("LastTileYSize", minTileSizeY);
  713.             newNode.SetAttributeValue("MinX", ulx);
  714.             newNode.SetAttributeValue("MaxX", lrx);
  715.             newNode.SetAttributeValue("MinY", lry);
  716.             newNode.SetAttributeValue("MaxY", uly);
  717.             newNode.SetAttributeValue("ScaleFactor", dec.ScaleX);
  718.             newNode.SetAttributeValue("TilesCountX", tilesCountX);
  719.             newNode.SetAttributeValue("TilesCountY", tilesCountY);
  720.  
  721.             doc.Element("SMPyramid").Add(newNode);
  722.  
  723.             doc.Save(xmlPath);
  724.  
  725.             fhInputTile.Dispose();
  726.         }
  727.  
  728.         public static IEnumerable<string> GetTiles(string xmlFileName, int level, Envelope viewportData, bool debug = false)
  729.         {
  730.             // load level data
  731.             XDocument doc = XDocument.Load(xmlFileName);
  732.             if (doc == null)
  733.             {
  734.                 throw new FileNotFoundException("File '" + xmlFileName + "' not found.");
  735.             }
  736.  
  737.             XElement levelData =
  738.                 doc.Element("SMPyramid")
  739.                     .Descendants()
  740.                     .First(element => int.Parse(element.Attribute("NR").Value) == level);
  741.  
  742.             XElement pyramid = doc.Element("SMPyramid");
  743.  
  744.             string pyramidPath = pyramid.Attribute("LayerPath").Value;
  745.             string layerName = pyramid.Attribute("LayerName").Value;
  746.             string extension = pyramid.Attribute("Extension").Value;
  747.  
  748.             double minX = double.Parse(levelData.Attribute("MinX").Value);
  749.             double maxX = double.Parse(levelData.Attribute("MaxX").Value);
  750.             double minY = double.Parse(levelData.Attribute("MinY").Value);
  751.             double maxY = double.Parse(levelData.Attribute("MaxY").Value);
  752.  
  753.             int tilesCountX = int.Parse(levelData.Attribute("TilesCountX").Value);
  754.             int tilesCountY = int.Parse(levelData.Attribute("TilesCountY").Value);
  755.  
  756.             double scaleFactor = double.Parse(levelData.Attribute("ScaleFactor").Value);
  757.             double rasterXSize = double.Parse(levelData.Attribute("RasterXSize").Value);
  758.             double rasterYSize = double.Parse(levelData.Attribute("RasterYSize").Value);
  759.             double lastTileXSize = double.Parse(levelData.Attribute("LastTileXSize").Value);
  760.             double lastTileYSize = double.Parse(levelData.Attribute("LastTileYSize").Value);
  761.  
  762.             Envelope[][] extents = new Envelope[tilesCountX][];
  763.  
  764.             double tileGeoWidth = rasterXSize * scaleFactor;
  765.             double tileGeoHeight = rasterYSize * scaleFactor;
  766.             for (int i = 0; i < extents.Length; i++)
  767.             {
  768.                 double offsetY = (double) i*tileGeoHeight;
  769.                 extents[i] = new Envelope[tilesCountY];
  770.                 for (int j = 0; j < extents[i].Length; j++)
  771.                 {
  772.                     // count extents for pieces
  773.                     double offsetX = (double) j*tileGeoWidth;
  774.                    
  775.                     extents[i][j] = new Envelope();
  776.                     extents[i][j].MinX = minX + offsetX;
  777.                     extents[i][j].MinY = maxY - tileGeoHeight - offsetY;
  778.                     extents[i][j].MaxX = minX + tileGeoWidth + offsetX;
  779.                     extents[i][j].MaxY = maxY - offsetY;
  780.  
  781.                     // last tiles are cut
  782.                     if (j == extents[i].Length - 1 && i == extents.Length - 1)
  783.                     {
  784.                         extents[i][j].MaxX -= ((rasterXSize - lastTileXSize) * scaleFactor);
  785.                         extents[i][j].MinY += ((rasterYSize - lastTileYSize) * scaleFactor);
  786.                     }
  787.                     else if (j == extents[i].Length - 1)
  788.                     {
  789.                         extents[i][j].MaxX -= ((rasterXSize - lastTileXSize) * scaleFactor);
  790.                     }
  791.                     else if (i == extents.Length - 1)
  792.                     {
  793.                         extents[i][j].MinY += ((rasterYSize - lastTileYSize) * scaleFactor);
  794.                     }
  795.                 }
  796.             }
  797.  
  798.             var intersectedIndices =
  799.                 from a in Enumerable.Range(0, tilesCountX)
  800.                 from b in Enumerable.Range(0, tilesCountY)
  801.                 where EnvIntersects(extents[a][b], viewportData)
  802.                 select new {a, b};
  803.  
  804.             foreach (var index in intersectedIndices)
  805.             {
  806.                 string[] pathElements =
  807.                     {
  808.                         pyramidPath,
  809.                         level.ToString(),
  810.                         index.a.ToString(),
  811.                         string.Format("{0}_{1}_{2}{3}", layerName, (index.a + 1).ToString().PadLeft(5, '0'), (index.b + 1).ToString().PadLeft(5, '0'), extension)
  812.                     };
  813.                 yield return Path.Combine(pathElements);
  814.             }
  815.  
  816.             if (debug)
  817.             {
  818.                 Console.WriteLine("GLOBAL_EXTENT = {0}, {1}, {2}, {3}", minX, maxX, minY, maxY);
  819.                 Console.WriteLine("TILES_COUNT = {0}/{1}", tilesCountX, tilesCountY);
  820.                 for (int i = 0; i < extents.LongLength; i++)
  821.                 {
  822.                     for (int j = 0; j < extents.Length; j++)
  823.                     {
  824.                         Console.WriteLine("EXTENTS = ({0}, {1}, {2}, {3})", extents[i][j].MinX, extents[i][j].MaxX, extents[i][j].MinY, extents[i][j].MaxY);
  825.                     }
  826.                 }
  827.                 Console.WriteLine("BOUNDS = ({0}, {1}, {2}, {3})", viewportData.MinX, viewportData.MaxX, viewportData.MinY, viewportData.MaxY);
  828.             }
  829.         }
  830.  
  831.         private static bool EnvIntersects(Envelope e1, Envelope e2)
  832.         {
  833.             return ((e1.MinX >= e2.MinX && e1.MinX <= e2.MaxX) || (e1.MinX < e2.MinX && e1.MaxX >= e2.MinX))
  834.                 && ((e1.MinY >= e2.MinY && e1.MinY <= e2.MaxY) || (e1.MinY < e2.MinY && e1.MaxY >= e2.MaxY));
  835.         }
  836.  
  837.         private DataSource BuildPyramidLevel(MosaicInfo levelMosaicInfo, TileInfo levelOutputTileInfo, int level,
  838.             Action<int> progressCallback, out int minTileSizeXOut, out int minTileSizeYOut, out string extension)
  839.         {
  840.             List<int> xRange = Enumerable.Range(1, levelOutputTileInfo.CountTilesX).ToList();
  841.             List<int> yRange = Enumerable.Range(1, levelOutputTileInfo.CountTilesY).ToList();
  842.  
  843.             DataSource ogrds = CreateTileIndex("TileResult_" + level, _tileIndexFieldName, _sourceSRS,
  844.                                                _tileIndexDriverTyp);
  845.  
  846.             int minTileSizeX = _tileWidth;
  847.             int minTileSizeY = _tileHeight;
  848.             string tileName = null;
  849.             int tileCounter = 0;
  850.             foreach (int yIndex in yRange)
  851.             {
  852.                 foreach (int xIndex in xRange)
  853.                 {
  854.                     int offsetY = (yIndex - 1)*levelOutputTileInfo.TileHeight;
  855.                     int offsetX = (xIndex - 1)*levelOutputTileInfo.TileWidth;
  856.  
  857.                     int width, height;
  858.                     if (yIndex == levelOutputTileInfo.CountTilesY)
  859.                     {
  860.                         height = levelOutputTileInfo.LastTileHeight;
  861.                     }
  862.                     else
  863.                     {
  864.                         height = levelOutputTileInfo.TileHeight;
  865.                     }
  866.  
  867.                     if (xIndex == levelOutputTileInfo.CountTilesX)
  868.                     {
  869.                         width = levelOutputTileInfo.LastTileWidth;
  870.                     }
  871.                     else
  872.                     {
  873.                         width = levelOutputTileInfo.TileWidth;
  874.                     }
  875.  
  876.                     tileName = GetTileName(levelMosaicInfo, levelOutputTileInfo, xIndex, yIndex, level);
  877.                     CreatePyramidTile(levelMosaicInfo, offsetX, offsetY, width, height, tileName, ogrds);
  878.  
  879.                     if (width < minTileSizeX)
  880.                     {
  881.                         minTileSizeX = width;
  882.                     }
  883.  
  884.                     if (height < minTileSizeY)
  885.                     {
  886.                         minTileSizeY = height;
  887.                     }
  888.  
  889.                     if (progressCallback != null)
  890.                     {
  891.                         if (yIndex == yRange.Last() && xIndex == xRange.Last())
  892.                         {
  893.                             progressCallback(100);
  894.                         }
  895.                         else
  896.                         {
  897.                             int progress =
  898.                                 (int) (((double) (++tileCounter)/(double) (xRange.Count*yRange.Count))*100d)%101;
  899.                             progressCallback(progress);
  900.                         }
  901.                     }
  902.                 }
  903.             }
  904.  
  905.             extension = Path.GetExtension(tileName);
  906.  
  907.             if (!string.IsNullOrEmpty(_tileIndexName))
  908.             {
  909.                 // TODO: check this
  910.                 string shapeName = GetTargetDir(level) + _tileIndexName;
  911.                 CopyTileIndexToDisk(ogrds, shapeName);
  912.             }
  913.  
  914.             if (!string.IsNullOrEmpty(_csvFileName))
  915.             {
  916.                 string csvName = GetTargetDir(level) + _csvFileName;
  917.                 CopyTileIndexToCSV(ogrds, csvName);
  918.             }
  919.  
  920.             minTileSizeXOut = minTileSizeX;
  921.             minTileSizeYOut = minTileSizeY;
  922.  
  923.             return ogrds;
  924.         }
  925.  
  926.         private void CreatePyramidTile(MosaicInfo levelMosaicInfo, int offsetX, int offsetY, int width, int height, string tileName, DataSource ogrds)
  927.         {
  928.             double sx = levelMosaicInfo.ScaleX * _scaleFactor;
  929.             double sy = levelMosaicInfo.ScaleY * _scaleFactor;
  930.  
  931.             AffineTransformDecorator dec =
  932.                 new AffineTransformDecorator(new[]
  933.                                                  {
  934.                                                      levelMosaicInfo.Ulx + offsetX*sx, sx, 0,
  935.                                                      levelMosaicInfo.Uly + offsetY*sy, 0, sy
  936.                                                  });
  937.  
  938.             Dataset sFh = levelMosaicInfo.GetDataSet((int) Math.Round(dec.Ulx),
  939.                                                      (int) Math.Round((dec.Uly + height*dec.ScaleY)),
  940.                                                      (int) Math.Round((dec.Ulx + width*dec.ScaleX)),
  941.                                                      (int) Math.Round(dec.Uly));
  942.  
  943.             if (sFh == null)
  944.             {
  945.                 return;
  946.             }
  947.  
  948.             if (ogrds != null)
  949.             {
  950.                 var points = dec.PointsFor(width, height);
  951.                 AddFeature(ogrds, tileName, points);
  952.             }
  953.  
  954.             DataType bt;
  955.             if (_bandType == DataType.GDT_Unknown)
  956.             {
  957.                 bt = levelMosaicInfo.BandType;
  958.             }
  959.             else
  960.             {
  961.                 bt = _bandType;
  962.             }
  963.  
  964.             double[] geotransform = { dec.Ulx, dec.ScaleX, 0, dec.Uly, 0, dec.ScaleY };
  965.             int bands = levelMosaicInfo.Bands;
  966.  
  967.             Dataset tFh;
  968.             if (_memDriver == null)
  969.             {
  970.                 tFh = _driver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
  971.             }
  972.             else
  973.             {
  974.                 tFh = _memDriver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
  975.             }
  976.  
  977.             if (tFh == null)
  978.             {
  979.                 throw new Exception("Creation failed, terminating gdal_tile.");
  980.             }
  981.  
  982.             tFh.SetGeoTransform(geotransform);
  983.             tFh.SetProjection(levelMosaicInfo.Projection);
  984.  
  985.             for (int band = 1; band < bands + 1; band++)
  986.             {
  987.                 Band tBand = tFh.GetRasterBand(band);
  988.                 tBand.Fill(255, 0); // otherwise it is filled with 0
  989.                 tBand.SetNoDataValue(255); // set metadata
  990.              
  991.                 if (levelMosaicInfo.Ct != null)
  992.                 {
  993.                     tBand.SetRasterColorTable(levelMosaicInfo.Ct);
  994.                 }
  995.                 tBand.SetRasterColorInterpretation(levelMosaicInfo.Ci[band - 1]);
  996.             }
  997.  
  998.             var res = Gdal.ReprojectImage(sFh, tFh, null, null, _resamplingMethod, 0, 0, null, null);
  999.  
  1000.             if (res != 0)
  1001.             {
  1002.                 throw new Exception("Reprojection failed for " + tileName + ", error " + res + ".");
  1003.             }
  1004.  
  1005.             levelMosaicInfo.CloseDataSet(ref sFh);
  1006.  
  1007.             if (_memDriver != null)
  1008.             {
  1009.                 Dataset ttFh = _driver.CreateCopy(tileName, tFh, 0, _createOptions.ToArray(), null, null);
  1010.                 ttFh.Dispose();
  1011.             }
  1012.  
  1013.             tFh.Dispose();
  1014.         }
  1015.  
  1016.         private DataSource TileImage(MosaicInfo mInfo, TileInfo ti)
  1017.         {
  1018.             _lastRowIndex = -1;
  1019.             DataSource ogrds = CreateTileIndex("TileResult_0", _tileIndexFieldName, _sourceSRS, _tileIndexDriverTyp);
  1020.  
  1021.             for (int yIndex = 1; yIndex < ti.CountTilesY + 1; yIndex++)
  1022.             {
  1023.                 for (int xIndex = 1; xIndex < ti.CountTilesX + 1; xIndex++)
  1024.                 {
  1025.                     int width, height;
  1026.                     int offsetY = (yIndex - 1) * ti.TileHeight;
  1027.                     int offsetX = (xIndex - 1) * ti.TileWidth;
  1028.                     if (yIndex == ti.CountTilesY)
  1029.                     {
  1030.                         height = ti.LastTileHeight;
  1031.                     }
  1032.                     else
  1033.                     {
  1034.                         height = ti.TileHeight;
  1035.                     }
  1036.  
  1037.                     if (xIndex == ti.CountTilesX)
  1038.                     {
  1039.                         width = ti.LastTileWidth;
  1040.                     }
  1041.                     else
  1042.                     {
  1043.                         width = ti.TileWidth;
  1044.                     }
  1045.  
  1046.                     string tileName = string.Empty;
  1047.                     if (UseDirForEachRow)
  1048.                     {
  1049.                         tileName = GetTileName(mInfo, ti, xIndex, yIndex, 0);
  1050.                     }
  1051.                     else
  1052.                     {
  1053.                         tileName = GetTileName(mInfo, ti, xIndex, yIndex);
  1054.                     }
  1055.  
  1056.                     CreateTile(ref mInfo, offsetX, offsetY, width, height, tileName, ref ogrds);
  1057.                 }
  1058.             }
  1059.  
  1060.             if (!string.IsNullOrEmpty(_tileIndexName))
  1061.             {
  1062.                 string shapeName = string.Empty;
  1063.                 if (UseDirForEachRow && !PyramidOnly)
  1064.                 {
  1065.                     shapeName = GetTargetDir(0) + _tileIndexName;
  1066.                 }
  1067.                 else
  1068.                 {
  1069.                     shapeName = GetTargetDir() + _tileIndexName;
  1070.                 }
  1071.  
  1072.                 CopyTileIndexToDisk(ogrds, shapeName);
  1073.             }
  1074.  
  1075.             if (!string.IsNullOrEmpty(_csvFileName))
  1076.             {
  1077.                 string csvName = string.Empty;
  1078.                 if (UseDirForEachRow && !PyramidOnly)
  1079.                 {
  1080.                     csvName = GetTargetDir(0) + _csvFileName;
  1081.                 }
  1082.                 else
  1083.                 {
  1084.                     csvName = GetTargetDir() + _csvFileName;
  1085.                 }
  1086.  
  1087.                 CopyTileIndexToCSV(ogrds, csvName);
  1088.             }
  1089.  
  1090.             return ogrds;
  1091.         }
  1092.  
  1093.         private void CopyTileIndexToDisk(DataSource ogrds, string fileName)
  1094.         {
  1095.             DataSource shapeDS = CreateTileIndex(fileName, _tileIndexFieldName, ogrds.GetLayerByIndex(0).GetSpatialRef(),
  1096.                                                  "ESRI Shapefile");
  1097.             ogrds.GetLayerByIndex(0).ResetReading();
  1098.  
  1099.             while (true)
  1100.             {
  1101.                 Feature feature = ogrds.GetLayerByIndex(0).GetNextFeature();
  1102.                 if (feature == null)
  1103.                 {
  1104.                     break;
  1105.                 }
  1106.  
  1107.                 Feature newFeature = feature.Clone();
  1108.                 string basename = Path.GetFileName(feature.GetFieldAsString(0));
  1109.                 if (UseDirForEachRow)
  1110.                 {
  1111.                     string fExt = Path.GetExtension(basename);
  1112.                     basename = fExt + "/" + basename;
  1113.                 }
  1114.  
  1115.                 newFeature.SetField(0, basename);
  1116.                 shapeDS.GetLayerByIndex(0).CreateFeature(newFeature);
  1117.             }
  1118.  
  1119.             CloseTileIndex(shapeDS);
  1120.  
  1121.             ogrds.Dispose();
  1122.         }
  1123.  
  1124.         private void CopyTileIndexToCSV(DataSource ogrds, string fileName)
  1125.         {
  1126.             using (File.Open(fileName, FileMode.Create))
  1127.             {
  1128.             }
  1129.  
  1130.             ogrds.GetLayerByIndex(0).ResetReading();
  1131.  
  1132.             while (true)
  1133.             {
  1134.                 Feature feature = ogrds.GetLayerByIndex(0).GetNextFeature();
  1135.                 if (feature == null)
  1136.                 {
  1137.                     break;
  1138.                 }
  1139.  
  1140.                 string basename = Path.GetFileName(feature.GetFieldAsString(0));
  1141.                 if (UseDirForEachRow)
  1142.                 {
  1143.                     string fExt = Path.GetExtension(basename);
  1144.                     basename = fExt + "/" + basename;
  1145.                 }
  1146.  
  1147.                 using (StreamWriter sw = new StreamWriter(fileName, true))
  1148.                 {
  1149.                     sw.Write(basename);
  1150.                     Geometry geom = feature.GetGeometryRef();
  1151.                     Envelope coords = new Envelope();
  1152.                     geom.GetEnvelope(coords);
  1153.  
  1154.                     // TODO: check this
  1155.                     sw.Write(_csvDelimiter + coords.MinX);
  1156.                     sw.Write(_csvDelimiter + coords.MaxX);
  1157.                     sw.Write(_csvDelimiter + coords.MinY);
  1158.                     sw.Write(_csvDelimiter + coords.MaxY);
  1159.                     sw.Write(Environment.NewLine);
  1160.                 }
  1161.             }
  1162.  
  1163.             ogrds.Dispose();
  1164.         }
  1165.  
  1166.         private void CloseTileIndex(DataSource ogrDataSource)
  1167.         {
  1168.             ogrDataSource.Dispose();
  1169.         }
  1170.  
  1171.         private void CreateTile(ref MosaicInfo mInfo, int offsetX, int offsetY, int width, int height, string tileName, ref DataSource ogrds)
  1172.         {
  1173.             DataType bt;
  1174.             if (_bandType == DataType.GDT_Unknown)
  1175.             {
  1176.                 bt = mInfo.BandType;
  1177.             }
  1178.             else
  1179.             {
  1180.                 bt = _bandType;
  1181.             }
  1182.  
  1183.             AffineTransformDecorator dec = new AffineTransformDecorator(new[]
  1184.                                                                             {
  1185.                                                                                 mInfo.Ulx,
  1186.                                                                                 mInfo.ScaleX,
  1187.                                                                                 0d,
  1188.                                                                                 mInfo.Uly,
  1189.                                                                                 0d,
  1190.                                                                                 mInfo.ScaleY
  1191.                                                                             });
  1192.  
  1193.             Dataset sFh = mInfo.GetDataSet(
  1194.                 (int) Math.Round((dec.Ulx + offsetX*dec.ScaleX)),
  1195.                 (int) Math.Round((dec.Uly + offsetY*dec.ScaleY + height*dec.ScaleY)),
  1196.                 (int) Math.Round((dec.Ulx + offsetX*dec.ScaleX + width*dec.ScaleX)),
  1197.                 (int) Math.Round((dec.Uly + offsetY*dec.ScaleY)));
  1198.  
  1199.             if (sFh == null)
  1200.             {
  1201.                 return;
  1202.             }
  1203.  
  1204.             var geotransform = new[]
  1205.                                    {
  1206.                                        dec.Ulx + offsetX * dec.ScaleX,
  1207.                                        dec.ScaleX,
  1208.                                        0d,
  1209.                                        dec.Uly + offsetY * dec.ScaleY,
  1210.                                        0d,
  1211.                                        dec.ScaleY
  1212.                                    };
  1213.  
  1214.             if (ogrds != null)
  1215.             {
  1216.                 AffineTransformDecorator dec2 = new AffineTransformDecorator(geotransform);
  1217.                 var points = dec2.PointsFor(width, height);
  1218.                 AddFeature(ogrds, tileName, points);
  1219.             }
  1220.  
  1221.             int bands = mInfo.Bands;
  1222.  
  1223.             Dataset tFh;
  1224.             if (_memDriver == null)
  1225.             {
  1226.                 tFh = _driver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
  1227.             }
  1228.             else
  1229.             {
  1230.                 tFh = _memDriver.Create(tileName, width, height, bands, bt, _createOptions.ToArray());
  1231.             }
  1232.  
  1233.             if (tFh == null)
  1234.             {
  1235.                 throw new Exception("Creation failed, terminating gdal_tile.");
  1236.             }
  1237.  
  1238.             tFh.SetGeoTransform(geotransform);
  1239.             if (_sourceSRS != null)
  1240.             {
  1241.                 string wkt = string.Empty;
  1242.                 _sourceSRS.ExportToWkt(out wkt);
  1243.                 tFh.SetProjection(wkt);
  1244.             }
  1245.  
  1246.             int readX = Math.Min(sFh.RasterXSize, width);
  1247.             int readY = Math.Min(sFh.RasterYSize, height);
  1248.  
  1249.             Band tBand, sBand;
  1250.             for (int band = 1; band < bands + 1; band++)
  1251.             {
  1252.                 sBand = sFh.GetRasterBand(band);
  1253.                 tBand = tFh.GetRasterBand(band);
  1254.  
  1255.                 if (mInfo.Ct != null)
  1256.                 {
  1257.                     tBand.SetRasterColorTable(mInfo.Ct);
  1258.                 }
  1259.  
  1260.                 // TODO: check this
  1261.                 byte[] data = new byte[readX * readY];
  1262.                 sBand.ReadRaster(0, 0, readX, readY, data, readX, readY, 0, 0);
  1263.                 tBand.WriteRaster(0, 0, readX, readY, data, readX, readY, 0, 0);
  1264.             }
  1265.            
  1266.             mInfo.CloseDataSet(ref sFh);
  1267.  
  1268.             Dataset ttFh;
  1269.             if (_memDriver != null)
  1270.             {
  1271.                 ttFh = _driver.CreateCopy(tileName, tFh, 0, _createOptions.ToArray(), null, null);
  1272.                 ttFh.Dispose();
  1273.             }
  1274.  
  1275.             tFh.Dispose();
  1276.         }
  1277.  
  1278.         private string GetTileName(MosaicInfo mInfo, TileInfo ti, int xIndex, int yIndex, int level = -1)
  1279.         {
  1280. //            int max = ti.CountTilesX;
  1281. //            if (ti.CountTilesY > max)
  1282. //            {
  1283. //                max = ti.CountTilesY;
  1284. //            }
  1285.  
  1286.             const int countDigits = 5; // max.ToString().Length;
  1287.             string fileName = Path.GetFileNameWithoutExtension(_layerFileName).Replace("@", "");
  1288.             string extension = Path.GetExtension(_layerFileName).Replace("@", "");
  1289.  
  1290.             string format = string.Empty;
  1291.             if (UseDirForEachRow)
  1292.             {
  1293.                 //format = Path.Combine(GetTargetDir(level) + yIndex,
  1294.                 //                             fileName + "_%0" + countDigits + "i_%0" + countDigits + "i");
  1295.                 format = Path.Combine(GetTargetDir(level) + yIndex,
  1296.                                       fileName + "_" + yIndex.ToString().PadLeft(countDigits, '0') + "_" +
  1297.                                       xIndex.ToString().PadLeft(countDigits, '0'));
  1298.  
  1299.                 // see if there was a switch in the row, if so then create new dir for row
  1300.                 if (_lastRowIndex < yIndex)
  1301.                 {
  1302.                     _lastRowIndex = yIndex;
  1303.                     if (!Directory.Exists(GetTargetDir(level) + yIndex))
  1304.                     {
  1305.                         Directory.CreateDirectory(GetTargetDir(level) + yIndex);
  1306.                     }
  1307.                 }
  1308.             }
  1309.             else
  1310.             {
  1311.                 //format = GetTargetDir(level) + fileName + "_%0" + countDigits + "i_%0" + countDigits + "i";
  1312.                 format = GetTargetDir(level) + fileName + "_" + yIndex.ToString().PadLeft(countDigits, '0') + "_" + xIndex.ToString().PadLeft(countDigits, '0');
  1313.             }
  1314.  
  1315.             // check for the extension that should be used
  1316.             if (string.IsNullOrEmpty(extension))
  1317.             {
  1318.                 format += extension;
  1319.             }
  1320.             else
  1321.             {
  1322.                 format += "." + _extension;
  1323.             }
  1324.  
  1325.             return format;
  1326.         }
  1327.  
  1328.         private string GetTargetDir(int level = -1)
  1329.         {
  1330.             string directory = Path.Combine(_targetDir, Path.GetFileNameWithoutExtension(_layerFileName));
  1331.             if (!directory.EndsWith("\\"))
  1332.             {
  1333.                 directory = directory + "\\";
  1334.             }
  1335.             return level == -1 ? directory : (Path.Combine(directory, level.ToString()) + "\\");
  1336.         }
  1337.  
  1338.         private DataSource GetTileIndexFromFiles(IEnumerable<string> inputTiles, string driverType, bool? force8To32Bit = null)
  1339.         {
  1340.             DataSource ogrTileIndexDS = CreateTileIndex("TileIndex", _tileIndexFieldName, null, driverType);
  1341.             foreach (string inputTile in inputTiles)
  1342.             {
  1343.                 Dataset fhInputTile;
  1344.  
  1345.                 if ((force8To32Bit != null) ? force8To32Bit.Value : _8To32Bit)
  1346.                 {
  1347.                     TiffConverter tiffConverter = new TiffConverter(new[] {"TiffConverter", inputTile});
  1348.                     fhInputTile = tiffConverter.Convert8To24Bit();
  1349.                 }
  1350.                 else
  1351.                 {
  1352.                     fhInputTile = Gdal.Open(inputTile, Access.GA_ReadOnly);
  1353.                 }
  1354.  
  1355.                 if (fhInputTile == null)
  1356.                 {
  1357.                     return null;
  1358.                 }
  1359.  
  1360.                 // TODO: check this
  1361.                 double[] argout = new double[6];
  1362.                 fhInputTile.GetGeoTransform(argout);
  1363.                 AffineTransformDecorator dec = new AffineTransformDecorator(argout);
  1364.                 List<Point> points = dec.PointsFor(fhInputTile.RasterXSize, fhInputTile.RasterYSize);
  1365.  
  1366.                 if (_rasterSize == -1 || _meters == -1)
  1367.                 {
  1368.                     _meters = dec.ScaleX * fhInputTile.RasterXSize;
  1369.                     _rasterSize = fhInputTile.RasterXSize;
  1370.                 }
  1371.  
  1372.                 AddFeature(ogrTileIndexDS, inputTile, points);
  1373.                 fhInputTile.Dispose();
  1374.             }
  1375.  
  1376.             if (_verbose)
  1377.             {
  1378.                 Debug.WriteLine("Finished");
  1379.             }
  1380.  
  1381.             return ogrTileIndexDS;
  1382.         }
  1383.  
  1384.         private void AddFeature(DataSource ogrDataSource, string location, List<Point> points)
  1385.         {
  1386.             Layer ogrLayer = ogrDataSource.GetLayerByIndex(0); // TODO: check this
  1387.             Feature ogrFeature = new Feature(ogrLayer.GetLayerDefn());
  1388.             if (ogrFeature == null)
  1389.             {
  1390.                 throw new Exception("Could not vreate the feature.");
  1391.             }
  1392.  
  1393.             ogrFeature.SetField(_tileIndexFieldName, location);
  1394.  
  1395.             string wkt = string.Format("POLYGON (({0} {1},{2} {3},{4} {5},{6} {7},{8} {9} ))",
  1396.                 points.ElementAt(0).X,
  1397.                 points.ElementAt(0).Y,
  1398.                 points.ElementAt(1).X,
  1399.                 points.ElementAt(1).Y,
  1400.                 points.ElementAt(2).X,
  1401.                 points.ElementAt(2).Y,
  1402.                 points.ElementAt(3).X,
  1403.                 points.ElementAt(3).Y,
  1404.                 points.ElementAt(0).X,
  1405.                 points.ElementAt(0).Y);
  1406.  
  1407.             var ogrGeometry = Ogr.CreateGeometryFromWkt(ref wkt, ogrLayer.GetSpatialRef());
  1408.             if (ogrGeometry == null)
  1409.             {
  1410.                 throw new Exception("Could not create geometry.");
  1411.             }
  1412.  
  1413.             ogrFeature.SetGeometryDirectly(ogrGeometry);
  1414.             ogrLayer.CreateFeature(ogrFeature);
  1415.             ogrFeature.Dispose(); // TODO: check this
  1416.         }
  1417.  
  1418.         private static DataSource CreateTileIndex(string dsName, string fieldName, SpatialReference srs, string driverName)
  1419.         {
  1420.             var ogrDriver = Ogr.GetDriverByName(driverName);
  1421.             if (ogrDriver == null)
  1422.             {
  1423.                 throw new Exception("ESRI Shapefile _driver not found.");
  1424.             }
  1425.  
  1426.             DataSource ogrDataSource = ogrDriver.Open(dsName, 1);
  1427.             if (ogrDataSource != null)
  1428.             {
  1429.                 // python:
  1430.                 // OGRDataSource.Destroy()
  1431.                 ogrDataSource.Dispose(); // TODO: check this
  1432.                 ogrDriver.DeleteDataSource(dsName);
  1433.             }
  1434.  
  1435.             ogrDataSource = ogrDriver.CreateDataSource(dsName, new string[] { }); // TODO: check this
  1436.             if (ogrDataSource == null)
  1437.             {
  1438.                 throw new Exception("Could not open datasource " + dsName);
  1439.             }
  1440.  
  1441.             Layer ogrLayer = ogrDataSource.CreateLayer("index", srs, wkbGeometryType.wkbPolygon, new string[] { }); // TODO: check this
  1442.             if (ogrLayer == null)
  1443.             {
  1444.                 throw new Exception("Could not create layer.");
  1445.             }
  1446.  
  1447.             FieldDefn ogrFieldDefn = new FieldDefn(fieldName, FieldType.OFTString);
  1448.             if (ogrFieldDefn == null)
  1449.             {
  1450.                 throw new Exception("Could not create FieldDefn for " + fieldName);
  1451.             }
  1452.  
  1453.             ogrFieldDefn.SetWidth(256);
  1454.             if (ogrLayer.CreateField(ogrFieldDefn, 1) != 0) // TODO: check this
  1455.             {
  1456.                 throw new Exception("Could not create field for " + fieldName);
  1457.             }
  1458.  
  1459.             return ogrDataSource;
  1460.         }
  1461.  
  1462.         private static void RemoveEmptyDirectories(string startLocation)
  1463.         {
  1464.             foreach (var directory in Directory.GetDirectories(startLocation))
  1465.             {
  1466.                 RemoveEmptyDirectories(directory);
  1467.                 if (Directory.GetFiles(directory).Length == 0 && Directory.GetDirectories(directory).Length == 0)
  1468.                 {
  1469.                     Directory.Delete(directory, false);
  1470.                 }
  1471.             }
  1472.         }
  1473.  
  1474.         private static void Usage()
  1475.         {
  1476.             Console.WriteLine("Usage: gdal_retile ");
  1477.             Console.WriteLine("        [-v] [-co NAME=VALUE]* [-of out_format]");
  1478.             Console.WriteLine("        [-ps pixelWidth pixelHeight]");
  1479.             Console.WriteLine("        [-ot  {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/");
  1480.             Console.WriteLine("               CInt16/CInt32/CFloat32/CFloat64}]");
  1481.             Console.WriteLine("        [-tileIndex _tileIndexName [-tileIndexField fieldName]]");
  1482.             Console.WriteLine("        [-csv fileName [-csvDelim delimiter]]");
  1483.             Console.WriteLine("        [-s_srs srs_def]  [-pyramidOnly] -levels numberoflevels");
  1484.             Console.WriteLine("        [-r {near/bilinear/cubic/cubicspline/lanczos}]");
  1485.             Console.WriteLine("        [-useDirForEachRow]");
  1486.             Console.WriteLine("        [-layerFileName layer_file_name]");
  1487.             Console.WriteLine("        [-targetDir tile_directory input_files]");
  1488.  
  1489.         }
  1490.  
  1491.         private void UsageFormat()
  1492.         {
  1493.             Console.WriteLine("Valid formats:");
  1494.             int count = Gdal.GetDriverCount();
  1495.             for (int i = 0; i < count; i++)
  1496.             {
  1497.                 _driver = Gdal.GetDriver(i);
  1498.                 Console.WriteLine(_driver.ShortName);
  1499.             }
  1500.         }
  1501.  
  1502.         // a class providing some usefull methods for affine Transformations
  1503.         internal class AffineTransformDecorator
  1504.         {
  1505.             private double[] _geotransform;
  1506.             private double _scaleX, _scaleY, _ulx, _uly;
  1507.  
  1508.             public double ScaleX
  1509.             {
  1510.                 get { return _scaleX; }
  1511.                 set { _scaleX = value; }
  1512.             }
  1513.  
  1514.             public double ScaleY
  1515.             {
  1516.                 get { return _scaleY; }
  1517.                 set { _scaleY = value; }
  1518.             }
  1519.  
  1520.             public double Ulx
  1521.             {
  1522.                 get { return _ulx; }
  1523.                 set { _ulx = value; }
  1524.             }
  1525.  
  1526.             public double Uly
  1527.             {
  1528.                 get { return _uly; }
  1529.                 set { _uly = value; }
  1530.             }
  1531.  
  1532.             public AffineTransformDecorator(double[] transform)
  1533.             {
  1534.                 _geotransform = transform;
  1535.                 _scaleX = _geotransform[1];
  1536.                 _scaleY = _geotransform[5];
  1537.  
  1538.                 if (_scaleY > 0)
  1539.                     _scaleY *= -1;
  1540.  
  1541.                 _ulx = _geotransform[0];
  1542.                 _uly = _geotransform[3];
  1543.             }
  1544.  
  1545.             public List<Point> PointsFor(int width, int height)
  1546.             {
  1547.                 List<Point> points = new List<Point>();
  1548.                 int w = (int)_scaleX * width;
  1549.                 int h = (int)_scaleY * height;
  1550.  
  1551.                 points.Add(new Point((int) Math.Round(_ulx), (int) Math.Round(_uly)));
  1552.                 points.Add(new Point((int) Math.Round((_ulx + w)), (int) Math.Round(_uly)));
  1553.                 points.Add(new Point((int) Math.Round(_ulx + w), (int) Math.Round(_uly + h)));
  1554.                 points.Add(new Point((int) Math.Round(_ulx), (int) Math.Round(_uly + h)));
  1555.  
  1556.                 return points;
  1557.             }
  1558.         }
  1559.  
  1560.         // a class holding information about a GDAL file or a GDAL fileset
  1561.         internal class MosaicInfo : IDisposable
  1562.         {
  1563.             private Driver _tempDriver;
  1564.             private string _fileName;
  1565.             private DataSetCache _cache;
  1566.             private DataSource _ogrTileIndexDS;
  1567.             private int _bands;
  1568.             private DataType _bandType;
  1569.             private string _projection;
  1570.             private double _scaleX, _scaleY;
  1571.             private ColorTable _ct;
  1572.             private ColorInterp[] _ci;
  1573.             private double _ulx, _uly, _lrx, _lry;
  1574.             private double _xSize, _ySize;
  1575.  
  1576.             public double XSize
  1577.             {
  1578.                 get { return _xSize; }
  1579.                 set { _xSize = value; }
  1580.             }
  1581.  
  1582.             public double YSize
  1583.             {
  1584.                 get { return _ySize; }
  1585.                 set { _ySize = value; }
  1586.             }
  1587.  
  1588.             public double ScaleX
  1589.             {
  1590.                 get { return _scaleX; }
  1591.                 set { _scaleX = value; }
  1592.             }
  1593.  
  1594.             public double ScaleY
  1595.             {
  1596.                 get { return _scaleY; }
  1597.                 set { _scaleY = value; }
  1598.             }
  1599.  
  1600.             public string Projection
  1601.             {
  1602.                 get { return _projection; }
  1603.                 set { _projection = value; }
  1604.             }
  1605.  
  1606.             public string FileName
  1607.             {
  1608.                 get { return _fileName; }
  1609.                 set { _fileName = value; }
  1610.             }
  1611.  
  1612.             public DataType BandType
  1613.             {
  1614.                 get { return _bandType; }
  1615.                 set { _bandType = value; }
  1616.             }
  1617.  
  1618.             public double Ulx
  1619.             {
  1620.                 get { return _ulx; }
  1621.                 set { _ulx = value; }
  1622.             }
  1623.  
  1624.             public double Uly
  1625.             {
  1626.                 get { return _uly; }
  1627.                 set { _uly = value; }
  1628.             }
  1629.  
  1630.             public double Lrx
  1631.             {
  1632.                 get { return _lrx; }
  1633.                 set { _lrx = value; }
  1634.             }
  1635.  
  1636.             public double Lry
  1637.             {
  1638.                 get { return _lry; }
  1639.                 set { _lry = value; }
  1640.             }
  1641.  
  1642.             public int Bands
  1643.             {
  1644.                 get { return _bands; }
  1645.                 set { _bands = value; }
  1646.             }
  1647.  
  1648.             public ColorTable Ct
  1649.             {
  1650.                 get { return _ct; }
  1651.                 set { _ct = value; }
  1652.             }
  1653.  
  1654.             public ColorInterp[] Ci
  1655.             {
  1656.                 get { return _ci; }
  1657.                 set { _ci = value; }
  1658.             }
  1659.  
  1660.             public MosaicInfo(string fileName, DataSource inputDS, bool convert8To32Bit)
  1661.             {
  1662.                 _tempDriver = Gdal.GetDriverByName("MEM");
  1663.                 _fileName = fileName;
  1664.                 _cache = new DataSetCache(convert8To32Bit);
  1665.                 _ogrTileIndexDS = inputDS;
  1666.  
  1667.                 // TODO: check this
  1668.                 // python: self.ogrTileIndexDS.GetLayer().ResetReading()
  1669.                 _ogrTileIndexDS.GetLayerByIndex(0).ResetReading();
  1670.                 Feature feature = _ogrTileIndexDS.GetLayerByIndex(0).GetNextFeature();
  1671.                 string imgLocation = feature.GetFieldAsString(0);
  1672.  
  1673.                 Dataset fhInputTile = _cache.Get(imgLocation);
  1674.  
  1675.                 _bands = fhInputTile.RasterCount;
  1676.                 _bandType = fhInputTile.GetRasterBand(1).DataType;
  1677.                 _projection = fhInputTile.GetProjection();
  1678.  
  1679.                 double[] argout = new double[6];
  1680.                 fhInputTile.GetGeoTransform(argout);
  1681.  
  1682.                 AffineTransformDecorator dec = new AffineTransformDecorator(argout);
  1683.  
  1684.                 _scaleX = dec.ScaleX;
  1685.                 _scaleY = dec.ScaleY;
  1686.  
  1687.                 ColorTable ct = fhInputTile.GetRasterBand(1).GetRasterColorTable();
  1688.                 if (ct != null)
  1689.                 {
  1690.                     _ct = ct.Clone();
  1691.                 }
  1692.                 else
  1693.                 {
  1694.                     _ct = null;
  1695.                 }
  1696.  
  1697.                 List<ColorInterp> ciList = new List<ColorInterp>();
  1698.                 for (int i = 0; i < _bands; i++)
  1699.                 {
  1700.                     ciList.Add(ColorInterp.GCI_Undefined);
  1701.                 }
  1702.  
  1703.                 if (_ci == null)
  1704.                 {
  1705.                     _ci = new ColorInterp[_bands];
  1706.                 }
  1707.  
  1708.                 for (int iBand = 0; iBand < _bands; iBand++)
  1709.                 {
  1710.                     _ci[iBand] = fhInputTile.GetRasterBand(iBand + 1).GetRasterColorInterpretation();
  1711.                 }
  1712.  
  1713.                 Envelope extent = new Envelope();
  1714.                 _ogrTileIndexDS.GetLayerByIndex(0).GetExtent(extent, 1);
  1715.  
  1716.                 _ulx = extent.MinX;
  1717.                 _uly = extent.MaxY;
  1718.                 _lrx = extent.MaxX;
  1719.                 _lry = extent.MinY;
  1720.  
  1721.                 _xSize = (int) Math.Round((Math.Round((_lrx - _ulx)/_scaleX)));
  1722.                 _ySize = Math.Abs((int) Math.Round((Math.Round((_uly - _lry)/_scaleY))));
  1723.  
  1724.                 //fhInputTile.Dispose();
  1725.  
  1726.                 //Console.WriteLine("[" + inputDS.GetName() + "] - " + _xSize + "x" + _ySize + ": " + _ulx + ", " + _lrx + ", " + _lry + ", " + _uly);
  1727.             }
  1728.  
  1729.             public Dataset GetDataSet(double minX, double minY, double maxX, double maxY)
  1730.             {
  1731.                 _ogrTileIndexDS.GetLayerByIndex(0).ResetReading();
  1732.                 _ogrTileIndexDS.GetLayerByIndex(0).SetSpatialFilterRect(minX, minY, maxX, maxY);
  1733.  
  1734.                 List<Feature> features = new List<Feature>();
  1735.                 Envelope envelope = null;
  1736.  
  1737.                 while (true)
  1738.                 {
  1739.                     Feature feature = _ogrTileIndexDS.GetLayerByIndex(0).GetNextFeature();
  1740.                     if (feature == null)
  1741.                     {
  1742.                         break;
  1743.                     }
  1744.                     features.Add(feature);
  1745.                     if (envelope == null)
  1746.                     {
  1747.                         // TODO: check this
  1748.                         // python:
  1749.                         // envelope=feature.GetGeometryRef().GetEnvelope()
  1750.                         envelope = new Envelope();
  1751.                         feature.GetGeometryRef().GetEnvelope(envelope);
  1752.                     }
  1753.                     else
  1754.                     {
  1755.                         Envelope featureEnv = new Envelope();
  1756.                         feature.GetGeometryRef().GetEnvelope(featureEnv);
  1757.  
  1758.                         envelope.MinX = Math.Min(featureEnv.MinX, envelope.MinX);
  1759.                         envelope.MinY = Math.Min(featureEnv.MinY, envelope.MinY);
  1760.                         envelope.MaxX = Math.Max(featureEnv.MaxX, envelope.MaxX);
  1761.                         envelope.MaxY = Math.Max(featureEnv.MaxY, envelope.MaxY);
  1762.                     }
  1763.                 }
  1764.  
  1765.                 if (envelope == null)
  1766.                 {
  1767.                     return null;
  1768.                 }
  1769.  
  1770.                 // enlarge to query rect if necessairy
  1771.                 envelope.MinX = Math.Min(minX, envelope.MinX);
  1772.                 envelope.MinY = Math.Min(minY, envelope.MinY);
  1773.                 envelope.MaxX = Math.Max(maxX, envelope.MaxX);
  1774.                 envelope.MaxY = Math.Max(maxY, envelope.MaxY);
  1775.  
  1776.                 _ogrTileIndexDS.GetLayerByIndex(0).SetSpatialFilter(null);
  1777.  
  1778.                 // merge tiles
  1779.                 int resultSizeX = (int) Math.Round((Math.Ceiling((double) (maxX - minX)/(double) _scaleX)));
  1780.                 int resultSizeY = (int) Math.Round((Math.Ceiling((double) (minY - maxY)/(double) _scaleY)));
  1781.  
  1782.                 Dataset resultDS = _tempDriver.Create("TEMP", resultSizeX, resultSizeY, _bands, _bandType, new string[] { });
  1783.                 resultDS.SetGeoTransform(new[] { minX, _scaleX, 0, maxY, 0, _scaleY });
  1784.  
  1785.                 foreach (Feature feature in features)
  1786.                 {
  1787.                     string featureName = feature.GetFieldAsString(0);
  1788.                     Dataset sourceDS = _cache.Get(featureName);
  1789.  
  1790.                     double[] argout = new double[6];
  1791.                     sourceDS.GetGeoTransform(argout);
  1792.  
  1793.                     AffineTransformDecorator dec = new AffineTransformDecorator(argout);
  1794.  
  1795.                     // calculate read and write offsets
  1796.                     int readOffsetX = (int) Math.Round((Math.Round((double) (minX - dec.Ulx)/(double) _scaleX)));
  1797.                     int readOffsetY = (int) Math.Round((Math.Round((double) (maxY - dec.Uly)/(double) _scaleY)));
  1798.  
  1799.                     int writeOffsetX = 0;
  1800.                     if (readOffsetX < 0)
  1801.                     {
  1802.                         writeOffsetX = readOffsetX * (-1);
  1803.                         readOffsetX = 0;
  1804.                     }
  1805.                     int writeOffsetY = 0;
  1806.                     if (readOffsetY < 0)
  1807.                     {
  1808.                         writeOffsetY = readOffsetY * (-1);
  1809.                         readOffsetY = 0;
  1810.                     }
  1811.  
  1812.                     // calculate read and write dimensions
  1813.                     int readX = Math.Min(Math.Min(resultSizeX, sourceDS.RasterXSize - readOffsetX), resultSizeX - writeOffsetX);
  1814.                     if (readX <= 0)
  1815.                         continue;
  1816.  
  1817.                     int readY = Math.Min(Math.Min(resultSizeY, sourceDS.RasterYSize - readOffsetY), resultSizeY - writeOffsetY);
  1818.                     if (readY <= 0)
  1819.                         continue;
  1820.  
  1821.                     for (int bandNr = 1; bandNr < _bands + 1; bandNr++)
  1822.                     {
  1823.                         Band sBand = sourceDS.GetRasterBand(bandNr);
  1824.                         Band tBand = resultDS.GetRasterBand(bandNr);
  1825.  
  1826.                         if (_ct != null)
  1827.                         {
  1828.                             tBand.SetRasterColorTable(_ct);
  1829.                         }
  1830.                         tBand.SetRasterColorInterpretation(_ci[bandNr - 1]);
  1831.  
  1832.                         byte[] buffer = new byte[readX * readY];
  1833.                         sBand.ReadRaster(readOffsetX, readOffsetY, readX, readY, buffer, readX, readY, 0, 0);
  1834.                         tBand.WriteRaster(writeOffsetX, writeOffsetY, readX, readY, buffer, readX, readY, 0, 0);
  1835.                     }
  1836.  
  1837. //                    sourceDS.Dispose();
  1838.                 }
  1839.  
  1840.  
  1841.                 return resultDS;
  1842.             }
  1843.  
  1844.             public void CloseDataSet(ref Dataset memDS)
  1845.             {
  1846.                 memDS.Dispose();
  1847.             }
  1848.  
  1849.             public void Dispose()
  1850.             {
  1851.                 _cache.Dispose();
  1852.                 _ogrTileIndexDS.Dispose();
  1853.             }
  1854.         }
  1855.  
  1856.         internal class DataSetCache : IDisposable
  1857.         {
  1858.             private int _cacheSize;
  1859.             private Queue<string> _queue;
  1860.             private Dictionary<string, Dataset> _dict;
  1861.             private bool _8To32Bit;
  1862.  
  1863.             public DataSetCache(bool convert8to32bit)
  1864.             {
  1865.                 _8To32Bit = convert8to32bit;
  1866.                 _cacheSize = 8;
  1867.                 _queue = new Queue<string>();
  1868.                 _dict = new Dictionary<string, Dataset>();
  1869.             }
  1870.  
  1871.             public Dataset Get(string name)
  1872.             {
  1873.                 if (_dict.ContainsKey(name))
  1874.                     return _dict[name];
  1875.  
  1876.                 Dataset result;
  1877.                 if (_8To32Bit)
  1878.                 {
  1879.                     TiffConverter tiffConverter = new TiffConverter(new[] { "TiffConverter", name });
  1880.                     result = tiffConverter.Convert8To24Bit();
  1881.                 }
  1882.                 else
  1883.                 {
  1884.                     result = Gdal.Open(name, Access.GA_ReadOnly);
  1885.                 }
  1886.  
  1887.                 if (result == null)
  1888.                 {
  1889.                     throw new Exception("Error opening: " + name);
  1890.                 }
  1891.  
  1892.                 if (_queue.Count() == _cacheSize)
  1893.                 {
  1894.                     string toRemove = _queue.Dequeue();
  1895.                     _dict.Remove(toRemove);
  1896.                 }
  1897.  
  1898.                 _queue.Enqueue(name);
  1899.                 _dict.Add(name, result);
  1900.  
  1901. //                result.Dispose();
  1902.  
  1903.                 return _dict[name];
  1904.             }
  1905.  
  1906.             public void Dispose()
  1907.             {
  1908.                 foreach (var keyValuePair in _dict)
  1909.                 {
  1910.                     keyValuePair.Value.Dispose();
  1911.                 }
  1912.                 _dict.Clear();
  1913.                 _queue.Clear();
  1914.                 _dict = null;
  1915.                 _queue = null;
  1916.                 GC.Collect();
  1917.             }
  1918.         }
  1919.  
  1920.         // a class holding info how to tile
  1921.         internal class TileInfo
  1922.         {
  1923.             private int _tileWidth, _tileHeight;
  1924.             private int _countTilesX, _countTilesY;
  1925.             private int _lastTileWidth, _lastTileHeight;
  1926.  
  1927.             public int CountTilesX
  1928.             {
  1929.                 get { return _countTilesX; }
  1930.                 set { _countTilesX = value; }
  1931.             }
  1932.  
  1933.             public int CountTilesY
  1934.             {
  1935.                 get { return _countTilesY; }
  1936.                 set { _countTilesY = value; }
  1937.             }
  1938.  
  1939.             public int TileWidth
  1940.             {
  1941.                 get { return _tileWidth; }
  1942.                 set { _tileWidth = value; }
  1943.             }
  1944.  
  1945.             public int TileHeight
  1946.             {
  1947.                 get { return _tileHeight; }
  1948.                 set { _tileHeight = value; }
  1949.             }
  1950.  
  1951.             public int LastTileWidth
  1952.             {
  1953.                 get { return _lastTileWidth; }
  1954.                 set { _lastTileWidth = value; }
  1955.             }
  1956.  
  1957.             public int LastTileHeight
  1958.             {
  1959.                 get { return _lastTileHeight; }
  1960.                 set { _lastTileHeight = value; }
  1961.             }
  1962.  
  1963.             public TileInfo(double xSize, double ySize, int tileWidth, int tileHeight)
  1964.             {
  1965.                 _tileWidth = tileWidth;
  1966.                 _tileHeight = tileHeight;
  1967.                 _countTilesX = (int) (xSize/(double) tileWidth);
  1968.                 _countTilesY = (int) (ySize/(double) tileHeight);
  1969.                 _lastTileWidth = (int) (xSize - (double) _countTilesX*(double) tileWidth);
  1970.                 _lastTileHeight = (int) (ySize - (double) _countTilesY*(double) tileHeight);
  1971.  
  1972.                 if (_lastTileWidth > 0)
  1973.                 {
  1974.                     _countTilesX = _countTilesX + 1;
  1975.                 }
  1976.                 else
  1977.                 {
  1978.                     _lastTileWidth = tileWidth;
  1979.                 }
  1980.  
  1981.                 if (_lastTileHeight > 0)
  1982.                 {
  1983.                     _countTilesY = _countTilesY + 1;
  1984.                 }
  1985.                 else
  1986.                 {
  1987.                     _lastTileHeight = tileHeight;
  1988.                 }
  1989.             }
  1990.  
  1991.             public void Report()
  1992.             {
  1993.                 Console.WriteLine("tileWidth       " + _tileWidth);
  1994.                 Console.WriteLine("tileHeight      " + _tileHeight);
  1995.                 Console.WriteLine("countTilesX:    " + _countTilesX);
  1996.                 Console.WriteLine("countTilesY:    " + _countTilesY);
  1997.                 Console.WriteLine("lastTileWidth:  " + _lastTileWidth);
  1998.                 Console.WriteLine("lastTileHeight: " + _lastTileHeight);
  1999.             }
  2000.         }
  2001.     }
  2002. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement