// Copyright 2007: Christian Graefe

// Copyright 2008: Dan Brecht and Joel Wilson

//

// This file is part of SharpMap.

// SharpMap is free software; you can redistribute it and/or modify

// it under the terms of the GNU Lesser General Public License as published by

// the Free Software Foundation; either version 2 of the License, or

// (at your option) any later version.

//

// SharpMap is distributed in the hope that it will be useful,

// but WITHOUT ANY WARRANTY; without even the implied warranty of

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

// GNU Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License

// along with SharpMap; if not, write to the Free Software

// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

using System;

using System.Buffers;

using System.Collections.Generic;

using System.Drawing;

using System.Drawing.Imaging;

using System.IO;

using Common.Logging;

using GeoAPI;

using GeoAPI.Geometries;

using OSGeo.GDAL;

using GeoAPI.CoordinateSystems;

using GeoAPI.CoordinateSystems.Transformations;

using NetTopologySuite.Index.Strtree;

using SharpMap.Base;

using SharpMap.CoordinateSystems;

using SharpMap.Data;

using SharpMap.Extensions.Data;

using Geometry = GeoAPI.Geometries.IGeometry;

using SharpMap.Rendering.Thematics;

using Point = System.Drawing.Point;

using Polygon = GeoAPI.Geometries.IPolygon;

namespace SharpMap.Layers

{

/// <summary>

/// Gdal raster image layer

/// </summary>

/// <remarks>

/// <example>

/// <code lang="C#">

/// myMap = new SharpMap.Map(new System.Drawing.Size(500,250);<br/>

/// SharpMap.Layers.GdalRasterLayer layGdal = new SharpMap.Layers.GdalRasterLayer("Blue Marble", @"C:\data\bluemarble.ecw");<br/>

/// myMap.Layers.Add(layGdal);<br/>

/// myMap.ZoomToExtents();<br/>

/// </code>

/// </example>

/// </remarks>

[Serializable]

public class GdalRasterLayer : Layer, ICanQueryLayer

{

private static readonly ILog _logger = LogManager.GetLogger(typeof (GdalRasterLayer));

static GdalRasterLayer()

{

GdalConfiguration.ConfigureGdal();

}

/// <summary>

/// Static Rtree that contains the extent and connection strings for gdal raster layer sub data sets.

/// </summary>

private STRtree<string> _subDataSets;

private IGeometryFactory _factory;

/// <summary>

/// Gets or sets a value indicating the geometry factory to use when creating geometries.

/// </summary>

/// <returns>A geometry factory</returns>

protected IGeometryFactory Factory

{

get

{

return _factory ?? (_factory = GeometryServiceProvider.Instance.CreateGeometryFactory(SRID));

}

set { _factory = value; }

}

/// <summary>

/// The extent of the layer

/// </summary>

protected Envelope _envelope;

/// <summary>

/// The inner dataset object

/// </summary>

protected Dataset _gdalDataset;

/// <summary>

/// The size of the raster image

/// </summary>

protected Size _imageSize;

// outer radius is feather between inner radius and rest of image

// outer radius is feather between inner radius and rest of image

private Point _stretchPoint;

/// <summary>

/// Flag indicating of geographic rotation is to be used.

/// </summary>

protected bool _useRotation;

#region accessors

/// <summary>

/// Gets the version of fwTools that was used to compile and test this GdalRasterLayer

/// </summary>

[Obsolete("FWTools no longer used", true)]

public static string FWToolsVersion

{

#pragma warning disable 612,618

get { return FwToolsHelper.FwToolsVersion; }

#pragma warning restore 612,618

}

/// <summary>

/// Gets or sets the filename of the raster file

/// </summary>

public string Filename { get; protected internal set; }

/// <summary>

/// Gets or sets a rectangle that is used to tile the image when rendering

/// </summary>

public Size TilingSize { get; set; }

/// <summary>

/// Gets or sets the bit depth of the raster file

/// </summary>

public int BitDepth { get; set; }

/// <summary>

/// Gets or set the projection of the raster file

/// </summary>

public string Projection { get; set; }

/// <summary>

/// Gets or sets to display IR Band

/// </summary>

public bool DisplayIR { get; set; }

/// <summary>

/// Gets or sets to display color InfraRed

/// </summary>

public bool DisplayCIR { get; set; }

/// <summary>

/// Gets or sets to display clip

/// </summary>

public bool ShowClip { get; set; }

/// <summary>

/// Gets or sets to display gamma

/// </summary>

public double Gamma { get; set; }

/// <summary>

/// Gets or sets to display gamma for Spot spot

/// </summary>

public double SpotGamma { get; set; }

/// <summary>

/// Gets or sets to display gamma for NonSpot

/// </summary>

public double NonSpotGamma { get; set; }

/// <summary>

/// Gets or sets to display red Gain

/// </summary>

public double[] Gain { get; set; }

/// <summary>

/// Gets or sets to display red Gain for Spot spot

/// </summary>

public double[] SpotGain { get; set; }

/// <summary>

/// Gets or sets to display red Gain for Spot spot

/// </summary>

public double[] NonSpotGain { get; set; }

/// <summary>

/// Gets or sets to display curve lut

/// </summary>

public List<int[]> CurveLut { get; set; }

/// <summary>

/// Correct Spot spot

/// </summary>

public bool HaveSpot { get; set; }

/// <summary>

/// Gets or sets to display curve lut for Spot spot

/// </summary>

public List<int[]> SpotCurveLut { get; set; }

/// <summary>

/// Gets or sets to display curve lut for NonSpot

/// </summary>

public List<int[]> NonSpotCurveLut { get; set; }

/// <summary>

/// Gets or sets the center point of the Spot spot

/// </summary>

public PointF SpotPoint { get; set; }

/// <summary>

/// Gets or sets the inner radius for the spot

/// </summary>

public double InnerSpotRadius { get; set; }

/// <summary>

/// Gets or sets the outer radius for the spot (feather zone)

/// </summary>

public double OuterSpotRadius { get; set; }

/// <summary>

/// Gets the true histogram

/// </summary>

public List<int[]> Histogram { get; private set; }

/// <summary>

/// Gets the quick histogram mean

/// </summary>

public double[] HistoMean

{

get

{

return null;// _histoMean;

}

}

/// <summary>

/// Gets the quick histogram brightness

/// </summary>

public double HistoBrightness

{

get

{

return double.NaN;

//_histoBrightness;

}

}

/// <summary>

/// Gets the quick histogram contrast

/// </summary>

public double HistoContrast

{

get

{

return double.NaN;

//_histoContrast;

}

}

/// <summary>

/// Gets the number of bands

/// </summary>

public int Bands { get; internal set; }

/// <summary>

/// Gets the GSD (Horizontal)

/// </summary>

public double GSD

{

get { return new GeoTransform(_gdalDataset).HorizontalPixelResolution; }

}

/// <summary>

/// Gets a value indicating if rotation information is to be used

/// </summary>

/// <returns><c>true</c> if rotation information is to be used.</returns>

public bool UseRotation

{

get { return _useRotation; }

set

{

if (value == _useRotation)

return;

_useRotation = value;

_envelope = GetExtent();

}

}

/// <summary>

/// Gets a value indicating the size of the raster data set.

/// </summary>

/// <returns>The size of the raster data set.</returns>

public Size Size

{

get { return _imageSize; }

}

/// <summary>

/// Gets or sets a value indicating if color correction is to be applied.

/// </summary>

public bool ColorCorrect { get; set; }

/// <summary>

/// Gets or sets a value indicating the portion of the image for which the <see cref="Histogram"/> is evaluated

/// </summary>

/// <returns>The portion of the image for which the <see cref="Histogram"/> is evaluated</returns>

public Rectangle HistoBounds { get; set; }

#pragma warning disable 1591

[Obsolete("Use CoordinateTransformation instead")]

public ICoordinateTransformation Transform

{

get { return CoordinateTransformation; }

protected set

{

CoordinateTransformation = value;

}

}

#pragma warning restore 1591

/// <summary>

/// Gets or sets a value indicating a color that is used as transparent color.

/// </summary>

/// <returns>A color that is interpreted as transparent.</returns>

public Color TransparentColor { get; set; }

/// <summary>

/// Gets or sets the minimum- and maximum pixel byte values of all raster bands.

/// <para/>Knowing these, you can scale the color range

/// </summary>

public Point StretchPoint

{

get

{

if (_stretchPoint.IsEmpty)

{

ComputeStretch();

}

return _stretchPoint;

}

set

{

_stretchPoint = value;

}

}

#endregion

/// <summary>

/// Creates an instance of this class

/// </summary>

/// <param name="layerName">The name of the layer</param>

protected GdalRasterLayer(string layerName)

{

SpotPoint = new PointF(0, 0);

Projection = "";

BitDepth = 8;

NonSpotGain = new double[] { 1, 1, 1, 1 };

SpotGain = new double[] { 1, 1, 1, 1 };

Gain = new double[] { 1, 1, 1, 1 };

NonSpotGamma = 1;

SpotGamma = 1;

Gamma = 1;

TransparentColor = Color.Empty;

ColorCorrect = true;

LayerName = layerName;

Gdal.AllRegister();

}

/// <summary>

/// initialize a Gdal based raster layer - from file

/// </summary>

/// <param name="strLayerName">Name of layer</param>

/// <param name="imageFilename">location of image</param>

public GdalRasterLayer(string strLayerName, string imageFilename) : this(strLayerName)

{

Filename = imageFilename;

OpenDataset(imageFilename);

}

/// <summary>

/// initialize a Gdal based raster layer - from buffer

/// </summary>

/// <param name="strLayerName">Name of layer</param>

/// <param name="imageBuffer">image buffer</param>

///

public GdalRasterLayer(string strLayerName, byte[] imageBuffer) : this(strLayerName)

{

// a layer data name

string memoryFileName = $"/vsimem/{strLayerName}Data";

Filename = memoryFileName;

OpenDataset(imageBuffer, memoryFileName);

}

/// <summary>

/// OpenDataset for buffer

/// </summary>

/// <param name="imageBuffer">image buffer</param>

/// <param name="memoryFileName">File Name in Memory</param>

protected void OpenDataset(byte[] imageBuffer, string memoryFileName)

{

try

{

Gdal.FileFromMemBuffer(memoryFileName, imageBuffer);

if (_logger.IsDebugEnabled)

_logger.Debug("Opening dataset " + memoryFileName);

_gdalDataset = Gdal.Open(memoryFileName, Access.GA_ReadOnly);

if (_logger.IsDebugEnabled)

_logger.Debug("Reading projection for " + memoryFileName);

// have gdal read the projection

Projection = _gdalDataset.GetProjectionRef();

if (!string.IsNullOrEmpty(Projection))

{

using (var p = new OSGeo.OSR.SpatialReference(null))

{

string wkt = Projection;

int res = p.ImportFromWkt(ref wkt);

int srid = 0;

if (res == 0)

{

if (Projection.StartsWith("PROJCS"))

int.TryParse(p.GetAuthorityCode("PROJCS"), out srid);

else if (Projection.StartsWith("GEOGCS"))

int.TryParse(p.GetAuthorityCode("GEOGCS"), out srid);

else

srid = p.AutoIdentifyEPSG();

}

SRID = srid;

}

}

if (_logger.IsDebugEnabled)

_logger.Debug("Read sizes" + memoryFileName);

_imageSize = new Size(_gdalDataset.RasterXSize, _gdalDataset.RasterYSize);

_envelope = GetExtent();

HistoBounds = new Rectangle((int)_envelope.MinX, (int)_envelope.MinY, (int)_envelope.Width,

(int)_envelope.Height);

//Bands = _gdalDataset.RasterCount;

if (_logger.IsDebugEnabled)

_logger.Debug("Dataset open, " + memoryFileName);

}

catch (Exception ex)

{

_gdalDataset = null;

throw new Exception("Couldn't load " + memoryFileName + "\n\n" + ex.Message + ex.InnerException);

}

finally

{

Gdal.Unlink(memoryFileName);

}

}

/// <summary>

/// OpenDataset for imageFile

/// </summary>

/// <param name="imageFilename">location of image</param>

protected void OpenDataset(string imageFilename)

{

try

{

if (_logger.IsDebugEnabled)

_logger.Debug("Opening dataset " + imageFilename);

_gdalDataset = Gdal.OpenShared(imageFilename, Access.GA_ReadOnly);

if (_logger.IsDebugEnabled)

_logger.Debug("Reading projection for " + imageFilename);

// have gdal read the projection

Projection = _gdalDataset.GetProjectionRef();

// no projection info found in the image...check for a prj

if (Projection == "" &&

File.Exists(imageFilename.Substring(0, imageFilename.LastIndexOf(".", StringComparison.CurrentCultureIgnoreCase)) + ".prj"))

{

Projection =

File.ReadAllText(imageFilename.Substring(0, imageFilename.LastIndexOf(".", StringComparison.CurrentCultureIgnoreCase)) + ".prj");

}

if (!String.IsNullOrEmpty(Projection))

{

using (var p = new OSGeo.OSR.SpatialReference(null))

{

var wkt = Projection;

var res = p.ImportFromWkt(ref wkt);

var srid = 0;

if (res == 0)

{

if (Projection.StartsWith("PROJCS"))

int.TryParse(p.GetAuthorityCode("PROJCS"), out srid);

else if (Projection.StartsWith("GEOGCS"))

int.TryParse(p.GetAuthorityCode("GEOGCS"), out srid);

else

srid = p.AutoIdentifyEPSG();

}

SRID = srid;

}

}

if (_logger.IsDebugEnabled)

_logger.Debug("Read sizes" + imageFilename);

_imageSize = new Size(_gdalDataset.RasterXSize, _gdalDataset.RasterYSize);

_envelope = GetExtent();

HistoBounds = new Rectangle((int) _envelope.MinX, (int) _envelope.MinY, (int) _envelope.Width,

(int) _envelope.Height);

//Bands = _gdalDataset.RasterCount;

if (_logger.IsDebugEnabled)

_logger.Debug("Dataset open, " + imageFilename);

}

catch (Exception ex)

{

_gdalDataset = null;

throw new Exception("Couldn't load " + imageFilename + "\n\n" + ex.Message + ex.InnerException);

}

}

/// <inheritdoc cref="DisposableObject.ReleaseManagedResources"/>

protected override void ReleaseManagedResources()

{

_factory = null;

if (_gdalDataset != null)

{

_gdalDataset.Dispose();

_gdalDataset = null;

}

base.ReleaseManagedResources();

}

/// <summary>

/// Returns the extent of the layer

/// </summary>

/// <returns>Bounding box corresponding to the extent of the features in the layer</returns>

public override Envelope Envelope

{

get { return _envelope.Copy(); }

}

/// <summary>

/// Renders the layer

/// </summary>

/// <param name="g">Graphics object reference</param>

/// <param name="map">Map which is rendered</param>

public override void Render(Graphics g, MapViewport map)

{

CheckDisposed();

ClearHistogram();

if (TilingSize.IsEmpty || (TilingSize.Width > map.Size.Width && TilingSize.Height > map.Size.Height))

{

if (_subDataSets != null)

{

foreach (string subDataSet in _subDataSets.Query(ToSource(map.Envelope)))

{

using (var ds = Gdal.Open(subDataSet, Access.GA_ReadOnly))

{

System.Diagnostics.Debug.WriteLine(ds.GetDescription());

GetPreview(ds, map.Size, g, map.Envelope, null, map);

}

}

}

else

GetPreview(_gdalDataset, map.Size, g, map.Envelope, null, map);

}

else

{

foreach (var tileEnvelope in Tile(map))

{

if (_subDataSets != null)

{

foreach (string subDataSet in _subDataSets.Query(ToSource(tileEnvelope)))

{

using (var ds = Gdal.Open(subDataSet, Access.GA_ReadOnly))

GetPreview(ds, map.Size, g, tileEnvelope, null, map);

}

}

else

GetPreview(_gdalDataset, map.Size, g, tileEnvelope, null, map);

}

}

//base.Render(g, map);

}

private void ClearHistogram()

{

var histogram = Histogram;

if (histogram == null)

{

histogram = new List<int[]>(Bands + 1);

for (int i = 0; i < Bands + 1; i++)

histogram.Add(new int[256]);

Histogram = histogram;

}

else

{

for (int i = 0; i < Bands + 1; i++)

Array.Clear(histogram[i], 0, 256);

}

}

private IEnumerable<Envelope> Tile(MapViewport map)

{

var lt = Point.Truncate(map.WorldToImage(map.Envelope.TopLeft()));

var rb = Point.Ceiling(map.WorldToImage(map.Envelope.BottomRight()));

var size = new Size(rb.X - lt.X, rb.Y - lt.Y);

var fullRect = new Rectangle(lt, size);

var overlapSize = new Size(2, 2);

for (int top = lt.Y; top < fullRect.Bottom; top+= TilingSize.Height )

{

for (int left = lt.X; left < fullRect.Right; left += TilingSize.Width)

{

var partialRect = new Rectangle(new Point(left - overlapSize.Width, top - overlapSize.Height),

Size.Add(TilingSize, overlapSize));

var res = new Envelope(map.ImageToWorld(partialRect.Location));

res.ExpandToInclude(map.ImageToWorld(new PointF(partialRect.Right, partialRect.Bottom)));

yield return res;

}

}

}

/// <summary>

/// Gets the spatial reference system of this layer

/// </summary>

/// <returns>A spatial reference system.</returns>

public ICoordinateSystem GetProjection()

{

try

{

return Session.Instance.CoordinateSystemServices.GetCoordinateSystem(SRID);

}

catch (Exception ee)

{

if (_logger.IsErrorEnabled)

_logger.Error("Error parsing projection", ee);

}

return null;

}

/// <summary>

/// Zoom to the native resolution

/// </summary>

/// <param name="map">The map object</param>

/// <returns>The zoom factor for a 1:1 Scale</returns>

public double GetOneToOne(Map map)

{

double dsWidth = _imageSize.Width;

double dsHeight = _imageSize.Height;

//double left, top, right, bottom;

double dblImgEnvW, dblImgEnvH, dblWindowGndW, dblWindowGndH, dblImginMapW, dblImginMapH;

var bbox = map.Envelope;

Size size = map.Size;

//// bounds of section of image to be displayed

//left = Math.Max(bbox.MinX, _envelope.MinX);

//top = Math.Min(bbox.MaxY, _envelope.MaxY);

//right = Math.Min(bbox.MaxX, _envelope.MaxX);

//bottom = Math.Max(bbox.MinY, _envelope.MinY);

// height and width of envelope of transformed image in ground space

dblImgEnvW = _envelope.Width;

dblImgEnvH = _envelope.Height;

// height and width of display window in ground space

dblWindowGndW = bbox.Width;

dblWindowGndH = bbox.Height;

// height and width of transformed image in pixel space

dblImginMapW = size.Width * (dblImgEnvW / dblWindowGndW);

dblImginMapH = size.Height * (dblImgEnvH / dblWindowGndH);

// image was not turned on its side

if ((dblImginMapW > dblImginMapH && dsWidth > dsHeight) ||

(dblImginMapW < dblImginMapH && dsWidth < dsHeight))

return map.Zoom * (dblImginMapW / dsWidth);

// image was turned on its side

else

return map.Zoom * (dblImginMapH / dsWidth);

}

/// <summary>

/// Calculates the Zoom value to nearest tiff internal resolution set

/// </summary>

/// <param name="map">A map</param>

/// <param name="bZoomIn">A flag indicating if zooming direction is magnify</param>

/// <returns>A zoom value</returns>

public double GetZoomNearestRSet(Map map, bool bZoomIn)

{

double dsWidth = _imageSize.Width;

double dsHeight = _imageSize.Height;

//double left, top, right, bottom;

double dblImgEnvW, dblImgEnvH, dblWindowGndW, dblWindowGndH, dblImginMapW, dblImginMapH;

double dblTempWidth;

Envelope bbox = map.Envelope;

Size size = map.Size;

//// bounds of section of image to be displayed

//left = Math.Max(bbox.MinX, _envelope.MinX);

//top = Math.Min(bbox.MaxY, _envelope.MaxY);

//right = Math.Min(bbox.MaxX, _envelope.MaxX);

//bottom = Math.Max(bbox.MinY, _envelope.MinY);

// height and width of envelope of transformed image in ground space

dblImgEnvW = _envelope.Width;

dblImgEnvH = _envelope.Height;

// height and width of display window in ground space

dblWindowGndW = bbox.Width;

dblWindowGndH = bbox.Height;

// height and width of transformed image in pixel space

dblImginMapW = size.Width * (dblImgEnvW / dblWindowGndW);

dblImginMapH = size.Height * (dblImgEnvH / dblWindowGndH);

// image was not turned on its side

if ((dblImginMapW > dblImginMapH && dsWidth > dsHeight) ||

(dblImginMapW < dblImginMapH && dsWidth < dsHeight))

dblTempWidth = dblImginMapW;

else

dblTempWidth = dblImginMapH;

// zoom level is within the r sets

if (dsWidth > dblTempWidth && (dsWidth / Math.Pow(2, 8)) < dblTempWidth)

{

if (bZoomIn)

{

for (int i = 0; i <= 8; i++)

{

if (dsWidth / Math.Pow(2, i) > dblTempWidth)

{

if (dsWidth / Math.Pow(2, i + 1) < dblTempWidth)

return map.Zoom * (dblTempWidth / (dsWidth / Math.Pow(2, i)));

}

}

}

else

{

for (int i = 8; i >= 0; i--)

{

if (dsWidth / Math.Pow(2, i) < dblTempWidth)

{

if (dsWidth / Math.Pow(2, i - 1) > dblTempWidth)

return map.Zoom * (dblTempWidth / (dsWidth / Math.Pow(2, i)));

}

}

}

}

return map.Zoom;

}

/// <summary>

/// Resets the <see cref="Histogram"/> bounds.

/// </summary>

public void ResetHistoRectangle()

{

HistoBounds = new Rectangle((int)_envelope.MinX, (int)_envelope.MinY,

(int)_envelope.Width, (int)_envelope.Height);

}

/// <summary>

/// Gets transform between raster's native projection and the map projection and sets it to <see cref="Layer.CoordinateTransformation"/>

/// </summary>

/// <param name="mapProjection">A map projection</param>

private void GetTransform(ICoordinateSystem mapProjection)

{

if (mapProjection == null || Projection == "")

{

CoordinateTransformation = null;

return;

}

// get our two projections

var srcCoord = GetProjection();

var tgtCoord = mapProjection;

// raster and map are in same projection, no need to transform

if (srcCoord.WKT == tgtCoord.WKT)

{

CoordinateTransformation = null;

return;

}

// create transform

CoordinateTransformation = Session.Instance.CoordinateSystemServices.CreateTransformation(srcCoord, tgtCoord);

ReverseCoordinateTransformation = Session.Instance.CoordinateSystemServices.CreateTransformation(tgtCoord, srcCoord);

}

private Envelope GetExtent()

{

var metadataDomainList = new HashSet<string>(_gdalDataset.GetMetadataDomainList());

if (metadataDomainList.Contains("SUBDATASETS"))

{

// Hack to clear number of bands

Bands = 0;

var envelope = new Envelope();

string[] subDataSets = _gdalDataset.GetMetadata("SUBDATASETS");

_subDataSets = new STRtree<string>();

for (int i = 0; i < subDataSets.Length; i++)

{

int namePos = subDataSets[i].IndexOf("_NAME=", StringComparison.Ordinal);

if (namePos < 0) continue;

string subDataSet = subDataSets[i].Substring(namePos + 6);

using (var ds = Gdal.Open(subDataSet, Access.GA_ReadOnly))

{

var tmp = GetExtent(ds);

_subDataSets.Insert(tmp, ds.GetDescription());

envelope.ExpandToInclude(tmp);

}

}

return envelope;

}

return GetExtent(_gdalDataset);

}

private Envelope GetExtent(Dataset dataSet)

{

if (dataSet == null)

return null;

// no rotation...use default transform

var geoTransform = new GeoTransform(dataSet);

if (geoTransform.IsIdentity)

geoTransform = new GeoTransform(-0.5, dataSet.RasterYSize + 0.5);

// image pixels

double dblW = dataSet.RasterXSize;

double dblH = dataSet.RasterYSize;

double left = geoTransform.EnvelopeLeft(dblW, dblH);

double right = geoTransform.EnvelopeRight(dblW, dblH);

double top = geoTransform.EnvelopeTop(dblW, dblH);

double bottom = geoTransform.EnvelopeBottom(dblW, dblH);

// Hack to set number of bands

if (Bands == 0)

Bands = dataSet.RasterCount;

return new Envelope(left, right, bottom, top);

}

/// <summary>

/// Gets the four corner coordinates of the raster data set.

/// </summary>

/// <returns>An array four corner coordinates of the raster data set</returns>

public Coordinate[] GetFourCorners()

{

var points = new Coordinate[4];

if (_gdalDataset != null)

{

var geoTrans = new double[6];

_gdalDataset.GetGeoTransform(geoTrans);

// no rotation...use default transform

if (!_useRotation && !HaveSpot || (DoublesAreEqual(geoTrans[0],0) && DoublesAreEqual(geoTrans[3],0)))

geoTrans = new[] { 999.5, 1, 0, 1000.5, 0, -1 };

points[0] = new Coordinate(geoTrans[0], geoTrans[3]);

points[1] = new Coordinate(geoTrans[0] + (geoTrans[1] * _imageSize.Width),

geoTrans[3] + (geoTrans[4] * _imageSize.Width));

points[2] = new Coordinate(geoTrans[0] + (geoTrans[1] * _imageSize.Width) + (geoTrans[2] * _imageSize.Height),

geoTrans[3] + (geoTrans[4] * _imageSize.Width) + (geoTrans[5] * _imageSize.Height));

points[3] = new Coordinate(geoTrans[0] + (geoTrans[2] * _imageSize.Height),

geoTrans[3] + (geoTrans[5] * _imageSize.Height));

// transform to map's projection

if (CoordinateTransformation != null)

{

for (var i = 0; i < 4; i++)

{

double[] dblPoint = CoordinateTransformation.MathTransform.Transform(new[] { points[i].X, points[i].Y });

points[i] = new Coordinate(dblPoint[0], dblPoint[1]);

}

}

}

return points;

}

/// <summary>

/// Gets a polygon outline of the raster data set boundary

/// </summary>

/// <returns>A polygon outline of the raster data set boundary</returns>

public Polygon GetFootprint()

{

var myRing = Factory.CreateLinearRing(GetFourCorners());

return Factory.CreatePolygon(myRing, null);

}

// applies map projection transfrom to get reprojected envelope

private void ApplyTransformToEnvelope()

{

_envelope = ToTarget(GetExtent(_gdalDataset));

HistoBounds = new Rectangle(

(int) Math.Floor(_envelope.MinX), (int) Math.Floor(_envelope.MinY),

(int) Math.Ceiling(_envelope.Width), (int) Math.Ceiling(_envelope.Width));

}

/// <summary>

/// Reprojects current raster to the given spatial reference system.

/// </summary>

/// <param name="cs">A spatial reference system.</param>

public void ReprojectToCoordinateSystem(ICoordinateSystem cs)

{

GetTransform(cs);

ApplyTransformToEnvelope();

}

// public method to set envelope and transform to new projection

/// <summary>

/// Method to set <see cref="Envelope"/> and <see cref="Layer.CoordinateTransformation"/> to the projection of the map

/// </summary>

/// <param name="map">The map</param>

public void ReprojectToMap(Map map)

{

ICoordinateSystem cs = null;

if (map.SRID > 0)

{

using (var p = new OSGeo.OSR.SpatialReference(null))

{

string wkt;

p.ImportFromEPSG(map.SRID);

p.ExportToWkt(out wkt);

cs = map.GetCoordinateSystem();

}

}

ReprojectToCoordinateSystem(cs);

}

/// <summary>

/// Renders the preview of the raster

/// </summary>

/// <param name="dataset">The raster data set</param>

/// <param name="size">Teh size</param>

/// <param name="g">A graphics object</param>

/// <param name="displayBbox">The bounding box of the display</param>

/// <param name="mapProjection">The spatial reference system of the map</param>

/// <param name="map">The map viewport</param>

protected virtual void GetPreview(Dataset dataset, Size size, Graphics g,

Envelope displayBbox, ICoordinateSystem mapProjection, MapViewport map)

{

//check if image is in bounding box

var dataSetBbox = ToTarget(GetExtent(dataset));

if (!displayBbox.Intersects(dataSetBbox))

return;

// Get the bounding box of the required data

dataSetBbox = ToSource(dataSetBbox.Intersection(displayBbox));

if (!NeedRotation(dataset))

{

GetNonRotatedPreview(dataset, size, g, displayBbox, mapProjection);

return;

}

var drawStart = DateTime.Now;

double totalReadDataTime = 0;

double totalTimeSetPixel = 0;