Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- using System;
- using System.Collections.Generic;
- using System.Diagnostics;
- using System.Drawing;
- using System.IO;
- using System.Linq;
- using System.Text.RegularExpressions;
- using System.Windows;
- using System.Windows.Media.Imaging;
- using System.Text;
- using System.Xml.Linq;
- using OSGeo.GDAL;
- using OSGeo.OGR;
- using OSGeo.OSR;
- using Driver = OSGeo.GDAL.Driver;
- using Point = System.Drawing.Point;
- using Size = System.Drawing.Size;
- 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 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 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;
- private string _layerFileName;
- private bool _pyramidOnly;
- private bool _useDirForEachRow;
- /// <summary>
- /// Creation options for output file. Multiple options can be specified.
- /// </summary>
- public enum CreationOptions
- {
- /// <summary>
- /// By default stripped TIFF files are created. This option can be used to force creation of tiled TIFF files.
- /// </summary>
- TILED,
- /// <summary>
- /// 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.
- /// </summary>
- GDAL_TIFF_INTERNAL_MASK,
- /// <summary>
- /// Set the compression to use: JPEG (JPEG files).
- /// </summary>
- COMPRESS_JPEG,
- /// <summary>
- /// Set the compression to use: LZW (PNG files).
- /// </summary>
- COMPRESS_LZW,
- /// <summary>
- /// Set the JPEG quality [1 - 100].
- /// </summary>
- JPEG_QUALITY,
- /// <summary>
- /// Set the JPEG quality [1 - 100].
- /// </summary>
- JPEG_QUALITY_OVERVIEW,
- /// <summary>
- /// 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).
- /// </summary>
- ALPHA
- }
- /// <summary>
- /// Output format, defaults to GeoTIFF (GTiff).
- /// </summary>
- public string OutputFormat
- {
- get { return _format; }
- set { _format = value; }
- }
- /// <summary>
- /// Force the output image bands to have a specific type. Use type names (ie. Byte, Int16,...).
- /// </summary>
- public string BandDataType
- {
- get { return Gdal.GetDataTypeName(_bandType); }
- set
- {
- _bandType = Gdal.GetDataTypeByName(value);
- if (_bandType == DataType.GDT_Unknown)
- {
- throw new Exception("Unknown GDAL data type: " + value);
- }
- }
- }
- /// <summary>
- /// 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).
- /// </summary>
- public string TargetDirectory
- {
- get { return _targetDir; }
- set
- {
- _targetDir = value;
- if (!Directory.Exists(_targetDir))
- {
- throw new Exception("TargetDir " + _targetDir + " does not exist.");
- }
- }
- }
- /// <summary>
- /// Pixel size to be used for the output file. If not specified, 512 x 512 is the default.
- /// </summary>
- public Size TileSize
- {
- get { return new Size(_tileWidth, _tileHeight); }
- set
- {
- _tileWidth = value.Width;
- _tileHeight = value.Height;
- }
- }
- /// <summary>
- /// Resampling algorithm, default is nearest-neighbour.
- /// </summary>
- public string ResamplingAlgorithm
- {
- get
- {
- switch (_resamplingMethod)
- {
- case ResampleAlg.GRA_NearestNeighbour:
- return "near";
- case ResampleAlg.GRA_Bilinear:
- return "bilinear";
- case ResampleAlg.GRA_Cubic:
- return "cubic";
- case ResampleAlg.GRA_CubicSpline:
- return "cubicspline";
- case (ResampleAlg)4:
- return "lanczos";
- default:
- throw new ArgumentOutOfRangeException();
- }
- }
- set
- {
- switch (value)
- {
- 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 ArgumentOutOfRangeException(value);
- }
- }
- }
- /// <summary>
- /// Number of pyramids levels to build (if not set, algorithm finishes, when a level consist of 1 tile only).
- /// </summary>
- public int Levels
- {
- get { return _levels; }
- set { _levels = value; }
- }
- /// <summary>
- /// 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).
- /// </summary>
- public string SpatialReference
- {
- get
- {
- string wkt = string.Empty;
- _sourceSRS.ExportToWkt(out wkt);
- return wkt;
- }
- set
- {
- string wkt = string.Empty;
- if (Osr.GetUserInputAsWKT(value, out wkt) != 0)
- {
- throw new Exception("Invalid -s_srs: " + value);
- }
- _sourceSRS = new SpatialReference(wkt);
- }
- }
- /// <summary>
- /// No retiling, build only the pyramids.
- /// </summary>
- public bool PyramidOnly
- {
- get { return _pyramidOnly; }
- set { _pyramidOnly = value; }
- }
- /// <summary>
- /// 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.
- /// </summary>
- public bool UseDirForEachRow
- {
- get { return _useDirForEachRow; }
- set { _useDirForEachRow = value; }
- }
- /// <summary>
- /// Convert 8 bit source image into 24 bit (recommended).
- /// </summary>
- public bool Convert8To24Bit
- {
- get { return _8To32Bit; }
- set { _8To32Bit = value; }
- }
- public string LayerFileName
- {
- get { return _layerFileName; }
- set { _layerFileName = value; }
- }
- /// <summary>
- /// Sets or removes chosen creation option.
- /// </summary>
- /// <param name="co">Chosen creation option.</param>
- /// <param name="add">If true, adds a new creation option. If false, removes it.</param>
- /// <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>
- public void SetCreationOption(CreationOptions co, bool add, string additionalValue = "")
- {
- string option = string.Empty;
- switch (co)
- {
- case CreationOptions.TILED:
- option = "TILED=YES";
- break;
- case CreationOptions.GDAL_TIFF_INTERNAL_MASK:
- option = "GDAL_TIFF_INTERNAL_MASK=TRUE";
- break;
- case CreationOptions.COMPRESS_JPEG:
- option = "COMPRESS=JPEG";
- break;
- case CreationOptions.COMPRESS_LZW:
- option = "COMPRESS=LZW";
- break;
- case CreationOptions.JPEG_QUALITY:
- option = "JPEG_QUALITY=" + additionalValue;
- break;
- case CreationOptions.JPEG_QUALITY_OVERVIEW:
- option = "JPEG_QUALITY_OVERVIEW=" + additionalValue;
- break;
- case CreationOptions.ALPHA:
- option = "ALPHA=YES";
- break;
- default:
- throw new ArgumentOutOfRangeException("Unknown creation option.");
- }
- if (add)
- _createOptions.Add(option);
- else
- _createOptions.Remove(option);
- }
- /// <summary>
- /// Sets input files for GdalRetile class.
- /// </summary>
- /// <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>
- /// <param name="isSourcePath">Specifies if the source is a path.</param>
- public void SetInputFiles(string source, bool isSourcePath = false)
- {
- string tmp = isSourcePath ? (new StreamReader(source)).ReadToEnd() : source;
- _names = tmp.Split(new[] { "\r\n", "\n" }, StringSplitOptions.RemoveEmptyEntries).ToList();
- }
- public GdalRetile()
- {
- Gdal.AllRegister();
- InitFields();
- }
- public GdalRetile(string[] args)
- {
- Gdal.AllRegister();
- InitFields();
- if (args.Any())
- {
- //string[] argv = Gdal.GeneralCmdLineProcessor(args, 0);
- string[] argv = args.Clone() as string[];
- 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));
- 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;
- case "-layerFileName":
- i++;
- _layerFileName = argv.ElementAt(i);
- break;
- default:
- if (arg.StartsWith("-") && !arg.Equals("--optfile"))
- {
- Usage();
- throw new Exception("Unrecognized command option: " + arg);
- }
- if (arg.Equals("--optfile"))
- {
- i++;
- using (StreamReader sr = new StreamReader(argv.ElementAt(i)))
- {
- _names = sr.ReadToEnd().Split(new[] {"\r\n", "\n"}, StringSplitOptions.RemoveEmptyEntries).ToList();
- }
- }
- else
- {
- _names.Add(arg);
- }
- break;
- }
- i++;
- }
- }
- else
- {
- throw new Exception("Command line has no parameters.");
- }
- }
- private void ValidateInputData()
- {
- 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 (string.IsNullOrEmpty(_layerFileName))
- {
- Usage();
- throw new Exception("Missing layer file name [-layerFileName].");
- }
- }
- public void MakePyramid(Action<int> progressCallback)
- {
- ValidateInputData();
- if (_format != "GTiff")
- {
- _memDriver = Gdal.GetDriverByName("MEM");
- }
- _driver = Gdal.GetDriverByName(_format);
- if (_driver == null)
- {
- UsageFormat();
- throw new Exception("Format _driver " + _format + " not found, pick a supported _driver.");
- }
- 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, false);
- if (tileIndexDS == null)
- {
- throw new Exception("Error building tile index.");
- }
- MosaicInfo mInfo = new MosaicInfo(_names[0], tileIndexDS, false);
- 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)
- {
- Action<int> callback;
- if (progressCallback != null)
- {
- callback = progressCallback;
- }
- else
- {
- callback = null;
- }
- BuildPyramid(dsCreatedTileIndex, _tileWidth, _tileHeight, callback);
- }
- }
- private void InitFields()
- {
- _verbose = false;
- _createOptions = new List<string>();
- _names = new List<string>();
- _tileWidth = 512;
- _tileHeight = 512;
- _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 = 10000;
- PyramidOnly = false;
- _lastRowIndex = -1;
- UseDirForEachRow = false;
- }
- private static double GetScale(int meters, int pixels)
- {
- return meters / pixels / MetersPerInch * DPI;
- }
- private void BuildPyramid(DataSource createdTileIndexDS, int tileWidth, int tileHeight, Action<int> progressCallback)
- {
- DataSource inputDS = createdTileIndexDS;
- for (int level = 0; level < _levels + 1; level++)
- {
- _lastRowIndex = -1;
- MosaicInfo levelMosaicInfo;
- switch (level)
- {
- case 0:
- _scaleFactor = 1;
- bool isJPEGCompression = (_createOptions.Any(option => option.Contains("COMPRESS=JPEG")))
- ? true
- : false;
- levelMosaicInfo = new MosaicInfo(_layerFileName, inputDS, isJPEGCompression); // original data for level 0 (for PNG)
- break;
- case 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;
- levelMosaicInfo = new MosaicInfo(_layerFileName, inputDS, _8To32Bit);
- break;
- default:
- _scaleFactor = 2;
- levelMosaicInfo = new MosaicInfo(_layerFileName, inputDS, _8To32Bit);
- break;
- }
- Action<int> callback;
- if (progressCallback != null)
- {
- callback = progressCallback;
- }
- else
- {
- callback = null;
- }
- TileInfo levelOutputTileInfo = new TileInfo(levelMosaicInfo.XSize/_scaleFactor,
- levelMosaicInfo.YSize/_scaleFactor, tileWidth, tileHeight);
- int minTileSizeX = 0, minTileSizeY = 0;
- string extension = string.Empty;
- inputDS = BuildPyramidLevel(levelMosaicInfo, levelOutputTileInfo, level, callback, out minTileSizeX, out minTileSizeY, out extension);
- StoreLevelData(level, inputDS, levelOutputTileInfo.CountTilesX, levelOutputTileInfo.CountTilesY, minTileSizeX, minTileSizeY, extension);
- levelMosaicInfo.Dispose();
- if (levelOutputTileInfo.CountTilesX * levelOutputTileInfo.CountTilesY == 1)
- {
- // remove empty dirs
- RemoveEmptyDirectories(Path.Combine(_targetDir, Path.GetFileNameWithoutExtension(_layerFileName)));
- return;
- }
- }
- }
- private void StoreLevelData(int nr, DataSource inputDS, int tilesCountX, int tilesCountY, int minTileSizeX, int minTileSizeY, string extension)
- {
- string xmlPath = Path.Combine(_targetDir, Path.GetFileName(_layerFileName));
- XDocument doc;
- if (nr == 0)
- {
- doc = new XDocument();
- XElement root = new XElement("SMPyramid");
- root.SetAttributeValue("LayerName", Path.GetFileNameWithoutExtension(_layerFileName));
- root.SetAttributeValue("LayerPath", Path.Combine(_targetDir, Path.GetFileNameWithoutExtension(_layerFileName)));
- root.SetAttributeValue("Extension", extension);
- doc.Add(root);
- }
- else
- {
- doc = XDocument.Load(xmlPath);
- }
- inputDS.GetLayerByIndex(0).ResetReading();
- Feature feature = inputDS.GetLayerByIndex(0).GetNextFeature();
- string imgLocation = feature.GetFieldAsString(0);
- Dataset fhInputTile = Gdal.Open(imgLocation, Access.GA_ReadOnly);
- double ulx = 0, uly = 0, lrx = 0, lry = 0;
- GdalHelper.GetGeocodedInfo(imgLocation, out ulx, out lrx, out lry, out uly);
- double[] argout = new double[6];
- fhInputTile.GetGeoTransform(argout);
- AffineTransformDecorator dec = new AffineTransformDecorator(argout);
- // Dataset lastTile = Gdal.Open(lastTileNameForLevel, Access.GA_ReadOnly);
- XElement newNode = new XElement("Level");
- newNode.SetAttributeValue("NR", nr);
- newNode.SetAttributeValue("RasterXSize", fhInputTile.RasterXSize);
- newNode.SetAttributeValue("RasterYSize", fhInputTile.RasterYSize);
- newNode.SetAttributeValue("LastTileXSize", minTileSizeX);
- newNode.SetAttributeValue("LastTileYSize", minTileSizeY);
- newNode.SetAttributeValue("MinX", ulx);
- newNode.SetAttributeValue("MaxX", lrx);
- newNode.SetAttributeValue("MinY", lry);
- newNode.SetAttributeValue("MaxY", uly);
- newNode.SetAttributeValue("ScaleFactor", dec.ScaleX);
- newNode.SetAttributeValue("TilesCountX", tilesCountX);
- newNode.SetAttributeValue("TilesCountY", tilesCountY);
- doc.Element("SMPyramid").Add(newNode);
- doc.Save(xmlPath);
- fhInputTile.Dispose();
- }
- public static IEnumerable<string> GetTiles(string xmlFileName, int level, Envelope viewportData, bool debug = false)
- {
- // load level data
- XDocument doc = XDocument.Load(xmlFileName);
- if (doc == null)
- {
- throw new FileNotFoundException("File '" + xmlFileName + "' not found.");
- }
- XElement levelData =
- doc.Element("SMPyramid")
- .Descendants()
- .First(element => int.Parse(element.Attribute("NR").Value) == level);
- XElement pyramid = doc.Element("SMPyramid");
- string pyramidPath = pyramid.Attribute("LayerPath").Value;
- string layerName = pyramid.Attribute("LayerName").Value;
- string extension = pyramid.Attribute("Extension").Value;
- double minX = double.Parse(levelData.Attribute("MinX").Value);
- double maxX = double.Parse(levelData.Attribute("MaxX").Value);
- double minY = double.Parse(levelData.Attribute("MinY").Value);
- double maxY = double.Parse(levelData.Attribute("MaxY").Value);
- int tilesCountX = int.Parse(levelData.Attribute("TilesCountX").Value);
- int tilesCountY = int.Parse(levelData.Attribute("TilesCountY").Value);
- double scaleFactor = double.Parse(levelData.Attribute("ScaleFactor").Value);
- double rasterXSize = double.Parse(levelData.Attribute("RasterXSize").Value);
- double rasterYSize = double.Parse(levelData.Attribute("RasterYSize").Value);
- double lastTileXSize = double.Parse(levelData.Attribute("LastTileXSize").Value);
- double lastTileYSize = double.Parse(levelData.Attribute("LastTileYSize").Value);
- Envelope[][] extents = new Envelope[tilesCountX][];
- double tileGeoWidth = rasterXSize * scaleFactor;
- double tileGeoHeight = rasterYSize * scaleFactor;
- for (int i = 0; i < extents.Length; i++)
- {
- double offsetY = (double) i*tileGeoHeight;
- extents[i] = new Envelope[tilesCountY];
- for (int j = 0; j < extents[i].Length; j++)
- {
- // count extents for pieces
- double offsetX = (double) j*tileGeoWidth;
- extents[i][j] = new Envelope();
- extents[i][j].MinX = minX + offsetX;
- extents[i][j].MinY = maxY - tileGeoHeight - offsetY;
- extents[i][j].MaxX = minX + tileGeoWidth + offsetX;
- extents[i][j].MaxY = maxY - offsetY;
- // last tiles are cut
- if (j == extents[i].Length - 1 && i == extents.Length - 1)
- {
- extents[i][j].MaxX -= ((rasterXSize - lastTileXSize) * scaleFactor);
- extents[i][j].MinY += ((rasterYSize - lastTileYSize) * scaleFactor);
- }
- else if (j == extents[i].Length - 1)
- {
- extents[i][j].MaxX -= ((rasterXSize - lastTileXSize) * scaleFactor);
- }
- else if (i == extents.Length - 1)
- {
- extents[i][j].MinY += ((rasterYSize - lastTileYSize) * scaleFactor);
- }
- }
- }
- var intersectedIndices =
- from a in Enumerable.Range(0, tilesCountX)
- from b in Enumerable.Range(0, tilesCountY)
- where EnvIntersects(extents[a][b], viewportData)
- select new {a, b};
- foreach (var index in intersectedIndices)
- {
- string[] pathElements =
- {
- pyramidPath,
- level.ToString(),
- index.a.ToString(),
- string.Format("{0}_{1}_{2}{3}", layerName, (index.a + 1).ToString().PadLeft(5, '0'), (index.b + 1).ToString().PadLeft(5, '0'), extension)
- };
- yield return Path.Combine(pathElements);
- }
- if (debug)
- {
- Console.WriteLine("GLOBAL_EXTENT = {0}, {1}, {2}, {3}", minX, maxX, minY, maxY);
- Console.WriteLine("TILES_COUNT = {0}/{1}", tilesCountX, tilesCountY);
- for (int i = 0; i < extents.LongLength; i++)
- {
- for (int j = 0; j < extents.Length; j++)
- {
- Console.WriteLine("EXTENTS = ({0}, {1}, {2}, {3})", extents[i][j].MinX, extents[i][j].MaxX, extents[i][j].MinY, extents[i][j].MaxY);
- }
- }
- Console.WriteLine("BOUNDS = ({0}, {1}, {2}, {3})", viewportData.MinX, viewportData.MaxX, viewportData.MinY, viewportData.MaxY);
- }
- }
- private static bool EnvIntersects(Envelope e1, Envelope e2)
- {
- return ((e1.MinX >= e2.MinX && e1.MinX <= e2.MaxX) || (e1.MinX < e2.MinX && e1.MaxX >= e2.MinX))
- && ((e1.MinY >= e2.MinY && e1.MinY <= e2.MaxY) || (e1.MinY < e2.MinY && e1.MaxY >= e2.MaxY));
- }
- private DataSource BuildPyramidLevel(MosaicInfo levelMosaicInfo, TileInfo levelOutputTileInfo, int level,
- Action<int> progressCallback, out int minTileSizeXOut, out int minTileSizeYOut, out string extension)
- {
- List<int> xRange = Enumerable.Range(1, levelOutputTileInfo.CountTilesX).ToList();
- List<int> yRange = Enumerable.Range(1, levelOutputTileInfo.CountTilesY).ToList();
- DataSource ogrds = CreateTileIndex("TileResult_" + level, _tileIndexFieldName, _sourceSRS,
- _tileIndexDriverTyp);
- int minTileSizeX = _tileWidth;
- int minTileSizeY = _tileHeight;
- string tileName = null;
- int tileCounter = 0;
- 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;
- }
- tileName = GetTileName(levelMosaicInfo, levelOutputTileInfo, xIndex, yIndex, level);
- CreatePyramidTile(levelMosaicInfo, offsetX, offsetY, width, height, tileName, ogrds);
- if (width < minTileSizeX)
- {
- minTileSizeX = width;
- }
- if (height < minTileSizeY)
- {
- minTileSizeY = height;
- }
- if (progressCallback != null)
- {
- if (yIndex == yRange.Last() && xIndex == xRange.Last())
- {
- progressCallback(100);
- }
- else
- {
- int progress =
- (int) (((double) (++tileCounter)/(double) (xRange.Count*yRange.Count))*100d)%101;
- progressCallback(progress);
- }
- }
- }
- }
- extension = Path.GetExtension(tileName);
- 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);
- }
- minTileSizeXOut = minTileSizeX;
- minTileSizeYOut = minTileSizeY;
- 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(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);
- tBand.Fill(255, 0); // otherwise it is filled with 0
- tBand.SetNoDataValue(255); // set metadata
- if (levelMosaicInfo.Ct != null)
- {
- tBand.SetRasterColorTable(levelMosaicInfo.Ct);
- }
- tBand.SetRasterColorInterpretation(levelMosaicInfo.Ci[band - 1]);
- }
- 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(ref 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(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(ref sFh);
- Dataset ttFh;
- if (_memDriver != null)
- {
- ttFh = _driver.CreateCopy(tileName, tFh, 0, _createOptions.ToArray(), null, null);
- ttFh.Dispose();
- }
- tFh.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;
- // }
- const int countDigits = 5; // max.ToString().Length;
- string fileName = Path.GetFileNameWithoutExtension(_layerFileName).Replace("@", "");
- string extension = Path.GetExtension(_layerFileName).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 = Path.Combine(_targetDir, Path.GetFileNameWithoutExtension(_layerFileName));
- if (!directory.EndsWith("\\"))
- {
- directory = directory + "\\";
- }
- return level == -1 ? directory : (Path.Combine(directory, level.ToString()) + "\\");
- }
- private DataSource GetTileIndexFromFiles(IEnumerable<string> inputTiles, string driverType, bool? force8To32Bit = null)
- {
- DataSource ogrTileIndexDS = CreateTileIndex("TileIndex", _tileIndexFieldName, null, driverType);
- foreach (string inputTile in inputTiles)
- {
- Dataset fhInputTile;
- if ((force8To32Bit != null) ? force8To32Bit.Value : _8To32Bit)
- {
- TiffConverter tiffConverter = new TiffConverter(new[] {"TiffConverter", inputTile});
- fhInputTile = tiffConverter.Convert8To24Bit();
- }
- else
- {
- 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(ogrTileIndexDS, inputTile, points);
- fhInputTile.Dispose();
- }
- if (_verbose)
- {
- Debug.WriteLine("Finished");
- }
- return ogrTileIndexDS;
- }
- private void AddFeature(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 RemoveEmptyDirectories(string startLocation)
- {
- foreach (var directory in Directory.GetDirectories(startLocation))
- {
- RemoveEmptyDirectories(directory);
- if (Directory.GetFiles(directory).Length == 0 && Directory.GetDirectories(directory).Length == 0)
- {
- Directory.Delete(directory, false);
- }
- }
- }
- 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(" [-layerFileName layer_file_name]");
- 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, bool convert8To32Bit)
- {
- _tempDriver = Gdal.GetDriverByName("MEM");
- _fileName = fileName;
- _cache = new DataSetCache(convert8To32Bit);
- _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();
- }
- 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))));
- //fhInputTile.Dispose();
- //Console.WriteLine("[" + inputDS.GetName() + "] - " + _xSize + "x" + _ySize + ": " + _ulx + ", " + _lrx + ", " + _lry + ", " + _uly);
- }
- 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]);
- 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);
- }
- // sourceDS.Dispose();
- }
- return resultDS;
- }
- public void CloseDataSet(ref 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;
- private bool _8To32Bit;
- public DataSetCache(bool convert8to32bit)
- {
- _8To32Bit = convert8to32bit;
- _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;
- if (_8To32Bit)
- {
- TiffConverter tiffConverter = new TiffConverter(new[] { "TiffConverter", name });
- result = tiffConverter.Convert8To24Bit();
- }
- else
- {
- 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.Add(name, result);
- // result.Dispose();
- return _dict[name];
- }
- public void Dispose()
- {
- foreach (var keyValuePair in _dict)
- {
- keyValuePair.Value.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) (xSize/(double) tileWidth);
- _countTilesY = (int) (ySize/(double) tileHeight);
- _lastTileWidth = (int) (xSize - (double) _countTilesX*(double) tileWidth);
- _lastTileHeight = (int) (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);
- }
- }
- }
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement