private void Identify(IEnumerable<ILayer> layers, Extent strict, Extent tolerant)
        {
            foreach (IMapLayer layer in layers)
            {
                IGroup grp = layer as IGroup;
                if (grp != null)
                {
                    Identify(grp, strict, tolerant);
                }
                else
                {
                    var gfl = layer as IMapFeatureLayer;
                    if (gfl != null)
                    {
                        if (gfl.DataSet.FeatureType == FeatureType.Polygon)
                        {
                            _frmFeatureIdentifier.Add(gfl, strict);
                        }
                        else
                        {
                            _frmFeatureIdentifier.Add(gfl, tolerant);
                        }
                    }

                    var rl = layer as IMapRasterLayer;
                    if (rl != null)
                    {
                        _frmFeatureIdentifier.Add(rl, strict);
                    }
                }
            }
        }

static double[] ToSequence(Extent extent)
{
    const int horizontal = 72;
    const int vertical = 36;
    var res = new double[horizontal * vertical * 2];

    var dx = extent.Width / (horizontal - 1);
    var dy = extent.Height / (vertical - 1);

    var minY = extent.MinY;
    var k = 0;
    for (var i = 0; i < vertical; i++)
    {
        var minX = extent.MinX;
        for (var j = 0; j < horizontal; j++)
        {
            res[k++] = minX;
            res[k++] = minY;
            minX += dx;
        }
        minY += dy;
    }

    return res;
}

        /// <summary>
        /// This function creates the HDR of Gridfile
        /// </summary>
        /// <param name="inExtents"> Extension of grid</param>
        /// <param name="cellSize">Size cell of the grid</param>
        /// <param name="projection">Projection (the same that shapefile)</param>
        /// <param name="noDataValue">No value definition</param>
        /// <param name="outGridPath">Path of the output</param>
        /// <param name="outGrid">Name of the output grid</param>
        public static void CreateGridFromExtents(
            Extent inExtents, double cellSize, ProjectionInfo projection, double noDataValue, string outGridPath, out IRaster outGrid)
        {

            double height = Math.Abs(inExtents.MaxY - inExtents.MinY);
            double width = Math.Abs(inExtents.MaxX - inExtents.MinX);
            int numberRows = Convert.ToInt32(Math.Ceiling(height / cellSize)) + 1;
            int numberCols = Convert.ToInt32(Math.Ceiling(width / cellSize)) + 1;

            // outGrid = Raster.CreateRaster(@outGridPath, null, demRaster.NumColumns, demRaster.NumRows, 1, demRaster.DataType, rasterOptions);
            outGrid = Raster.CreateRaster(@outGridPath, null, numberCols, numberRows, 1, typeof(float), new String[] { });
            outGrid.NoDataValue = noDataValue;
            outGrid.Projection = projection;
            outGrid.CellHeight = cellSize;
            outGrid.CellWidth = cellSize;
            //if (inExtents.MinX < 0)
            //    outGrid.Xllcenter = inExtents.MinX + (cellSize / 2.0);
            //else
            outGrid.Xllcenter = inExtents.MinX;// -(cellSize / 2.0);

            //if (inExtents.MinY < 0)
            //    outGrid.Yllcenter = inExtents.MinY + (cellSize / 2.0);
            //else
            outGrid.Yllcenter = inExtents.MinY;// -(cellSize / 2.0);
        }

 /// <summary>
 /// Creates a new raster bounds that is georeferenced to the specified envelope.
 /// </summary>
 /// <param name="numRows">The number of rows</param>
 /// <param name="numColumns">The number of columns</param>
 /// <param name="bounds">The bounding envelope</param>
 public RasterBounds(int numRows, int numColumns, Extent bounds)
 {
     _affine = new double[6];
     _numRows = numRows;
     _numColumns = numColumns;
     Extent = bounds;
 }

 /// <summary>
 /// The constructor for the ToolDialog
 /// </summary>
 /// <param name="tool">The ITool to create the dialog box for</param>
 /// <param name="dataSets">The list of available DataSets available</param>
 /// <param name="mapExtent">Creates a new instance of the tool dialog with map extent.</param>
 public ToolDialog(ITool tool, List<DataSetArray> dataSets, Extent mapExtent)
 {
     // Required by the designer
     InitializeComponent();
     DataSets = dataSets;
     _extent = mapExtent;
     Initialize(tool);
 }

 /// <summary>
 /// Converts a single geographic envelope into an equivalent Rectangle
 /// as it would be drawn on the screen.
 /// </summary>
 /// <param name="env">The geographic IEnvelope</param>
 /// <returns>A Rectangle</returns>
 public static Rectangle ProjToPixel(this IBasicMap map, Extent env)
 {
     var mapFrame = map.MapFrame as MapFrame;
     if (mapFrame != null) 
         return mapFrame.ProjToPixel(env);
     else
         return Rectangle.Empty;
 }

 /// <summary>
 /// Creates a new instance of an Extent Param with the specified name
 /// and the specified projection as the default projection that will
 /// appear if no changes are made.
 /// </summary>
 /// <param name="name"></param>
 /// <param name="defaultExtent"></param>
 public ExtentParam(string name, Extent defaultExtent)
 {
     Name = name;
     Value = defaultExtent;
     ParamVisible = ShowParamInModel.No;
     ParamType = "DotSpatial Extent Param";
     DefaultSpecified = true;
 }

static Extent Reproject(Extent extent, ProjectionInfo source, ProjectionInfo target, int depth = 0)
{
    var xy = ToSequence(extent);
    DotSpatial.Projections.Reproject.ReprojectPoints(xy, null, source, target, 0, xy.Length / 2);
    var res = ToExtent(xy);

    return res;
}

 public void SetAreaRectangle(Extent extent, ProjectionInfo rectangleProjection)
 {
     var xMin = extent.MinX;
     var yMin = extent.MinY;
     var xMax = extent.MaxX;
     var yMax = extent.MaxY;
     var box = new Box(xMin, xMax, yMin, yMax);
     SetAreaRectangle(box, rectangleProjection);
 }

        public LiDARDataSet(string filename)
        {
            //here read the maxX, maxY from the las file header
            //and change the Extent property
            LasReader Reader = new LasReader(filename);
            Reader.initialize();
            ulong PointNum = Reader.getNumPoints();
            int ArrayLength = (int)PointNum;
            //Range_double range = Reader.getBounds();
            //Reader.getBounds
            //Bounds_double ll = Reader.getBounds();
            Schema schema = Reader.getSchema();
            PointBuffer data = new PointBuffer(schema, (uint)PointNum);

            // get the dimensions (fields) of the point record for the X, Y, and Z values
            int offsetX = schema.getDimensionIndex(DimensionId.Id.X_i32);
            int offsetY = schema.getDimensionIndex(DimensionId.Id.Y_i32);
            int offsetZ = schema.getDimensionIndex(DimensionId.Id.Z_i32);
            Dimension dimensionX = schema.getDimension((uint)offsetX);
            Dimension dimensionY = schema.getDimension((uint)offsetY);
            Dimension dimensionZ = schema.getDimension((uint)offsetZ);

            // make the iterator to read from the file
            StageSequentialIterator iter = Reader.createSequentialIterator();
            uint numRead = iter.read(data);

            int xraw = data.getField_Int32(0, offsetX);
            int yraw = data.getField_Int32(0, offsetY);

            // LAS stores the data in scaled integer form: undo the scaling
            // so we can see the real values as doubles
            double MinX, MaxX, MinY, MaxY;
            MinX = MaxX = dimensionX.applyScaling_Int32(xraw);
            MinY = MaxY = dimensionY.applyScaling_Int32(yraw);

            for (int i = 1; i < ArrayLength; i++)
            {
                xraw = data.getField_Int32((uint)i, offsetX);
                yraw = data.getField_Int32((uint)i, offsetY);

                // LAS stores the data in scaled integer form: undo the scaling
                // so we can see the real values as doubles
                double x = dimensionX.applyScaling_Int32(xraw);
                double y = dimensionY.applyScaling_Int32(yraw);

                if (x < MinX) MinX = x;
                if (x > MaxX) MaxX = x;
                if (y < MinY) MinY = y;
                if (y > MaxY) MaxY = y;
            }
            setupExtent = new Extent(MinX, MinY, MaxX, MaxY);
        }

        /// <summary>
        /// The geographic envelope gives the region that the image should be created for.
        /// The window gives the corresponding pixel dimensions for the image, so that
        /// images matching the resolution of the screen can be used.
        /// </summary>
        /// <param name="envelope">
        /// The geographic extents to retrieve data for
        /// </param>
        /// <param name="window">
        /// The rectangle that defines the size of the drawing area in pixels
        /// </param>
        /// <returns>
        /// A bitmap captured from the main image
        /// </returns>
        public Bitmap GetBitmap(Extent envelope, Rectangle window)
        {
            if (window.Width == 0 || window.Height == 0)
                return null;

            if (Bounds == null || Bounds.Extent == null || Bounds.Extent.IsEmpty())
                return null;

            // Gets the scaling factor for converting from geographic to pixel coordinates
            double dx = (window.Width / envelope.Width);
            double dy = (window.Height / envelope.Height);

            double[] a = Bounds.AffineCoefficients;

            // gets the affine scaling factors.
            float m11 = Convert.ToSingle(a[1] * dx);
            float m22 = Convert.ToSingle(a[5] * -dy);
            float m21 = Convert.ToSingle(a[2] * dx);
            float m12 = Convert.ToSingle(a[4] * -dy);
            float l = (float)(a[0] - .5 * (a[1] + a[2])); // Left of top left pixel
            float t = (float)(a[3] - .5 * (a[4] + a[5])); // top of top left pixel
            float xShift = (float)((l - envelope.MinX) * dx);
            float yShift = (float)((envelope.MaxY - t) * dy);

            Bitmap tempResult = null;
            Bitmap result = null;
            Graphics g = null;
            try
            {
                tempResult = new Bitmap(window.Width, window.Height);
                g = Graphics.FromImage(tempResult);
                g.Transform = new Matrix(m11, m12, m21, m22, xShift, yShift);
                g.PixelOffsetMode = PixelOffsetMode.Half;
                if (m11 > 1 || m22 > 1)
                    g.InterpolationMode = InterpolationMode.NearestNeighbor;
                if (!g.VisibleClipBounds.IsEmpty)
                    g.DrawImage(_myImage, new PointF(0, 0));
                result = tempResult;
                tempResult = null;
            }
            catch (OverflowException) { } //Raised by g.DrawImage if the new images extent is to small
            finally
            {
                if (tempResult != null) tempResult.Dispose();
                if (g != null) g.Dispose();
            }
            return result;
        }

        /// <summary>
        /// This method simulates loading the array of points from the LAS file.
        /// Right now it is generating random points that are within the 
        /// view extent.
        /// </summary>
        /// <param name="boundingBox">the view extent</param>
        /// <returns>array of the points in [x y x y ... order]</returns>
        public double[] GetPointArray(Extent boundingBox)
        {
            double[] pointArray = new double[1000];
            Random rnd = new Random();
            double xMin = boundingBox.MinX;
            double yMin = boundingBox.MinY;

            for (int i = 0; i < 1000; i++)
            {
                double randomX = xMin + rnd.NextDouble() * boundingBox.Width;
                double randomY = yMin + rnd.NextDouble() * boundingBox.Height;
                pointArray[i] = randomX;
                i = i + 1;
                pointArray[i] = randomY;
            }
            return pointArray;
        }

        public SearchResult GetSeriesCatalogInPolygon(IList<IFeature> polygons, string[] keywords, double tileWidth, double tileHeight,
            DateTime startDate, DateTime endDate, WebServiceNode[] serviceIDs, BusinessObjects.Models.IProgressHandler bgWorker)
        {
            if (polygons == null) throw new ArgumentNullException("polygons");
            if (bgWorker == null) throw new ArgumentNullException("bgWorker");
            if (polygons.Count == 0)
            {
                throw new ArgumentException("The number of polygons must be greater than zero.");
            }

            if (keywords == null || keywords.Length == 0)
            {
                keywords = new[] { String.Empty };
            }

            var fullSeriesList = new List<BusinessObjects.Models.SeriesDataCartModel.SeriesDataCart>();
            for (int index = 0; index < polygons.Count; index++)
            {
                if (polygons.Count > 1)
                {
                    bgWorker.ReportMessage(string.Format("Processing polygons: {0} of {1}", index + 1, polygons.Count));
                }

                bgWorker.CheckForCancel();
                var polygon = polygons[index];
                var extentBox = new Extent(polygon.Envelope);
                var seriesForPolygon = GetSeriesListForExtent(extentBox, keywords, tileWidth, tileHeight, startDate,
                                                              endDate,
                                                              serviceIDs, bgWorker,
                                                              item => polygon.Intersects(new Coordinate(item.Longitude, item.Latitude)));
                fullSeriesList.AddRange(seriesForPolygon);
            }

            SearchResult resultFs = null;
            if (fullSeriesList.Count > 0)
            {
                bgWorker.ReportMessage("Calculating Points...");
                resultFs = SearchHelper.ToFeatureSetsByDataSource(fullSeriesList);
            }

            bgWorker.CheckForCancel();
            var message = string.Format("{0} Series found.", totalSeriesCount);
            bgWorker.ReportProgress(100, "Search Finished. " + message);
            return resultFs;
        }

        /// <summary>
        /// Gets the bitmap for the specified geographic envelope scaled to fit on a bitmap of the specified size in pixels.
        /// </summary>
        /// <param name="envelope"></param>
        /// <param name="pixelSize"></param>
        /// <returns></returns>
        public virtual Bitmap GetBitmap(Extent envelope, Size pixelSize)
        {
            Bitmap result = new Bitmap(pixelSize.Width, pixelSize.Height);
            Graphics g = Graphics.FromImage(result);
            foreach (ImageData image in _images)
            {
                Extent bounds = envelope.Intersection(image.Extent);

                Size ps = new Size((int)(pixelSize.Width * bounds.Width / envelope.Width),
                                   (int)(pixelSize.Height * bounds.Height / envelope.Height));
                int x = pixelSize.Width * (int)((bounds.X - envelope.X) / envelope.Width);
                int y = pixelSize.Height * (int)((envelope.Y - bounds.Y) / envelope.Height);
                if (ps.Width > 0 && ps.Height > 0)
                {
                    Bitmap tile = image.GetBitmap(bounds, ps);
                    g.DrawImageUnscaled(tile, x, y);
                }
            }
            return result;
        }

        /// <summary>
        /// Divides the search bounding box into several 'tiles' to prevent
        /// </summary>
        /// <param name="bigBoundingBox">the original bounding box</param>
        /// <param name="tileWidth">The tile width in decimal degrees</param>
        /// <param name="tileHeight">The tile height (south-north) in decimal degrees</param>
        /// <returns></returns>
        public static List<Extent> CreateTiles(Extent bigBoundingBox, double tileWidth, double tileHeight)
        {
            var tiles = new List<Extent>();
            double fullWidth = Math.Abs(bigBoundingBox.MaxX - bigBoundingBox.MinX);
            double fullHeight = Math.Abs(bigBoundingBox.MaxY - bigBoundingBox.MinY);

            if (fullWidth < tileWidth || fullHeight < tileHeight)
            {
                tiles.Add(bigBoundingBox);
                return tiles;
            }

            double yll = bigBoundingBox.MinY; //y-coordinate of the tile's lower left corner
            var numColumns = (int)(Math.Ceiling(fullWidth / tileWidth));
            var numRows = (int)(Math.Ceiling(fullHeight / tileHeight));
            var lastTileWidth = fullWidth - ((numColumns - 1) * tileWidth);
            var lastTileHeight = fullHeight - ((numRows - 1) * tileHeight);
            int r;

            for (r = 0; r < numRows; r++)
            {
                double xll = bigBoundingBox.MinX; //x-coordinate of the tile's lower left corner

                if (r == numRows - 1)
                {
                    tileHeight = lastTileHeight;
                }

                int c;
                for (c = 0; c < numColumns; c++)
                {
                    var newTile = c == (numColumns - 1) ? new Extent(xll, yll, xll + lastTileWidth, yll + tileHeight) :
                                                          new Extent(xll, yll, xll + tileWidth, yll + tileHeight);
                    tiles.Add(newTile);
                    xll = xll + tileWidth;
                }
                yll = yll + tileHeight;
            }
            return tiles;
        }

        public void GetNoDataCellCountTest()
        {
            string path = System.IO.Path.GetDirectoryName(System.Reflection.Assembly.GetExecutingAssembly().Location) + @"\..\..\..\Data\GetNoDataCellCountTest.BGD";

            const double xllcorner = 3267132.224761;
            const double yllcorner = 5326939.203029;
            const int ncols = 512;
            const int nrows = 128;
            const int frequencyOfNoValue = 5;

            const double cellsize = 500;
            double x2 = xllcorner + (cellsize * ncols);
            double y2 = yllcorner + (cellsize * nrows);
            Extent myExtent = new Extent(xllcorner, yllcorner, x2, y2);
            Raster target;
            target = Raster.Create(path, String.Empty, ncols, nrows, 1, typeof(double), new[] { String.Empty }) as Raster;
            target.Bounds = new RasterBounds(nrows, ncols, myExtent);
            target.NoDataValue = -9999;
            int mRow = target.Bounds.NumRows;
            int mCol = target.Bounds.NumColumns;

            for (int row = 0; row < mRow; row++)
            {
                for (int col = 0; col < mCol; col++)
                {
                    if (row % frequencyOfNoValue == 0)
                        target.Value[row, col] = -9999d;
                    else
                        target.Value[row, col] = 2d;
                }
            }
            target.Save();

            long expected = (nrows / frequencyOfNoValue) * ncols + ncols;
            long actual;
            actual = target.GetNoDataCellCount();
            Assert.AreEqual(expected, actual);

            System.IO.File.Delete(path);
        }

        public KriggingSimulator(Map map, IFeatureSet selectedDataset, Extent extent,
            Clip clip, string field, string layerName, string pixel)
        {
            this._field = field;
            this._map = map;
            this._selectedDataset = selectedDataset;
            this._layerName = layerName;
            this._extent = extent;
            this._clip = clip;
            this._pixel = pixel;

            this._prefixPath = this._field + "-" + Guid.NewGuid().ToString();
            this._folder = Path.Combine(AppDomain.CurrentDomain.BaseDirectory, "TempFiles");

            zedGraphControl1 = new ZedGraphControl();
            theoreticalModelControl1 = new TheoreticalModelControl();
            uxSearchNeighborhoodControl1 = new SearchNeighborhoodControl();
            uxOutputRaster1 = new OutputRaster();
            zedGraphControl2 = new ZedGraphControl();

            kriggingResult = new KriggingResult();
        }

        public void RasterMath_OutputHasSameBounds()
        {
            // Prepare input raster
            const double xllcorner = 3267132.224761;
            const double yllcorner = 5326939.203029;
            const int ncols = 39;
            const int nrows = 57;
            const double cellsize = 500;
            const double x2 = xllcorner + (cellsize * ncols);
            const double y2 = yllcorner + (cellsize * nrows);
            var myExtent = new Extent(xllcorner, yllcorner, x2, y2);
            var source = new Raster<int>(nrows, ncols)
            {
                Bounds = new RasterBounds(nrows, ncols, myExtent),
                NoDataValue = -9999
            };
            var mRow = source.Bounds.NumRows;
            var mCol = source.Bounds.NumColumns;

            var i = 0;
            for (var row = 0; row < mRow; row++)
            {
                for (var col = 0; col < mCol; col++)
                {
                    source.Value[row, col] = i++;
                }
            }

            var target = new RasterMultiply();
            IRaster outRaster = new Raster {Filename = FileTools.GetTempFileName(".bgd")};
            target.Execute(source, source, outRaster, new MockProgressHandler());
            outRaster = Raster.Open(outRaster.Filename);
            File.Delete(outRaster.Filename);

            Assert.AreEqual(source.NumColumns, outRaster.NumColumns);
            Assert.AreEqual(source.NumRows, outRaster.NumRows);
            Assert.AreEqual(source.Bounds.Extent, outRaster.Bounds.Extent);
        }

        /// <summary>
        /// The geographic envelope gives the region that the image should be created for.
        /// The window gives the corresponding pixel dimensions for the image, so that
        /// images matching the resolution of the screen can be used.
        /// </summary>
        /// <param name="envelope">
        /// The geographic extents to retrieve data for
        /// </param>
        /// <param name="window">
        /// The rectangle that defines the size of the drawing area in pixels
        /// </param>
        /// <returns>
        /// A bitmap captured from the main image
        /// </returns>
        public Bitmap GetBitmap(Extent envelope, Rectangle window)
        {
            if (window.Width == 0 || window.Height == 0)
            {
                return null;
            }

            Bitmap result = new Bitmap(window.Width, window.Height);
            Graphics g = Graphics.FromImage(result);

            //

            // Gets the scaling factor for converting from geographic to pixel coordinates
            double dx = (window.Width / envelope.Width);
            double dy = (window.Height / envelope.Height);

            double[] a = Bounds.AffineCoefficients;

            // gets the affine scaling factors.
            float m11 = Convert.ToSingle(a[1] * dx);
            float m22 = Convert.ToSingle(a[5] * -dy);
            float m21 = Convert.ToSingle(a[2] * dx);
            float m12 = Convert.ToSingle(a[4] * -dy);
            float l = (float)(a[0] - .5 * (a[1] + a[2])); // Left of top left pixel
            float t = (float)(a[3] - .5 * (a[4] + a[5])); // top of top left pixel
            float xShift = (float)((l - envelope.MinX) * dx);
            float yShift = (float)((envelope.MaxY - t) * dy);
            g.Transform = new Matrix(m11, m12, m21, m22, xShift, yShift);
            g.PixelOffsetMode = PixelOffsetMode.Half;
            if (m11 > 1 || m22 > 1)
            {
                g.InterpolationMode = InterpolationMode.NearestNeighbor;
            }

            g.DrawImage(_myImage, new PointF(0, 0));
            g.Dispose();
            return result;
        }

        /// <summary>
        /// Generates a new extent from the shape header.  This will infer the whether the ExtentMZ, ExtentM
        /// or Extent class is the best implementation.  Casting is required to access the higher
        /// values from the Extent return type.
        /// </summary>
        /// <returns>Extent, which can be Extent, ExtentM, or ExtentMZ</returns>
        public Extent ToExtent()
        {
            if (ShapeType == ShapeType.MultiPointZ ||
                ShapeType == ShapeType.PointZ ||
                ShapeType == ShapeType.PolygonZ ||
                ShapeType == ShapeType.PolyLineZ)
            {
                return new ExtentMZ(_xMin, _yMin, _mMin, _zMin, _xMax, _yMax, _mMax, _zMax);
            }
            if (ShapeType == ShapeType.MultiPointM ||
                ShapeType == ShapeType.PointM ||
                ShapeType == ShapeType.PolygonM ||
                ShapeType == ShapeType.PolyLineM)
            {
                return new ExtentM(_xMin, _yMin, _mMin, _xMax, _yMax, _mMax);
            }
            Extent ext = new Extent(_xMin, _yMin, _xMax, _yMax);

            return ext;
        }