Querying Raster Images

[David11]

Problem - We need to be able to return a raster cell value for a given latitude and longitude??

Background - We are currently evaluating ThinkGeo Mapsuite to use in a Land-use modelling environment, whereby we are running water-balance simulations (using CUDA) at defined gridpoints using a series of vector and raster layers as inputs to the models. We need to be able to query the raster layers to return cell values at the predefined grid point locations. Obviously this is easy to do for vector layers (although its a bit slow in the interregation), but we can't see how to do this for the raster layers. Currently our raster data is in the Erdas Image format using the FDO libraries, although we are likely to use other raster formats as well.

Possible Solutions - Here are the options:

    
1. Use the current mapsuite tools to do this?? Is this possible?? I can't see that it is.

    
2. Use the FDO or GDAL libraries to query the data - we will be looking into this next.

    
3. Modify the mapsuite tools to achieve this - again - is this possible?

Coming from a modelling environment, this is an extremely important functionality to have in a GIS tool... I had just assumed that the ThinkGeo mapsuite would be able to do this, and was quite shocked to find that it doesn't seem to have the capability. Does anyone have any suggestions or advice?

Regards, 

David McClymont

Principal Software Engineer

DHM Environmental Software Engineering Pty Ltd

[David]

David, 
  
   Thank you for your very precise e-mail.  I believe that we can do this without any changes.  I will send you a proof of concept shortly.  If this works then what I suggest is that we add a nice shortcut API to the RasterSource and RasterLayer base classes so this will be available to all of the inherited RasterLayers.  Our approach to this will be to create a very small rectangle around the point and request a 1x1 bitmap from the RasterSource.  After that we will get the cell value from the bitmap and pass it along.  I believe this will work.  When we implement the helper method we the Template Method Pattern as there are a number of raster layers we can easily get the cell value without using the trick above.  I will let you know our results. 
  
 David

[David11]

Thanks for your very prompt reply David. That would be awesome if you could show that this can be done. I'm sure there would be other modellers wanting to use MapSuite who would appreciate this functionality.

cheers

David

[David]

 David,

 
  Below is the code that should help you get the pixel value. If this lines up with your expectation we can integrate it into the daily builds so it will be tightly couples with the product.  We look forward to your suggestions and comments.
 
David
 



using System.Drawing;
using ThinkGeo.MapSuite.Core;

namespace ThinkGeo.MapSuite.Core
{
    static class RasterSourceHelper
    {
        public static GetPixelValueResult GetPixelValue(RasterSource rasterSource, PointShape pointShape)
        {
            GetPixelValueResult getPixelValueResult = new GetPixelValueResult();

            RectangleShape boudingBox = rasterSource.GetBoundingBox();

            // Here we times 1.0000000000001 to make sure the width/height is a little larger than a cell
            double pixelWidth = boudingBox.Width / rasterSource.GetImageWidth() * 1.0000000000001;
            double pixelHeight = boudingBox.Height / rasterSource.GetImageHeight() * 1.0000000000001;

            double minX = boudingBox.LowerLeftPoint.X;
            double minY = boudingBox.LowerLeftPoint.Y;
            double maxX = boudingBox.UpperRightPoint.X;
            double maxY = boudingBox.UpperRightPoint.Y;

            if (pointShape.X < minX || pointShape.X > maxX || pointShape.Y < minY || pointShape.Y > maxY)
            {
                // The input pointshape is not within the extent of rasterSource.
                getPixelValueResult.IsValid = false;
                return getPixelValueResult;
            }

            double smallRectangleShapeMinX = pointShape.X - pixelWidth;
            double smallRectangleShapeMaxY = pointShape.Y + pixelHeight;
            double smallRectangleShapeMaxX = pointShape.X;
            double smallRectangleShapeMinY = pointShape.Y;

            if (pointShape.X - pixelWidth < minX) // left
            {
                smallRectangleShapeMinX = pointShape.X;
                smallRectangleShapeMaxX = pointShape.X + pixelWidth;
            }
            else if (pointShape.X + pixelWidth > maxX)// right
            {
                smallRectangleShapeMinX = pointShape.X - pixelWidth;
                smallRectangleShapeMaxX = pointShape.X;
            }

            if (pointShape.Y + pixelHeight > maxY) // top
            {
                smallRectangleShapeMinY = pointShape.Y - pixelHeight;
                smallRectangleShapeMaxY = pointShape.Y;
            }
            else if (pointShape.Y - pixelHeight < minY) // bottom
            {
                smallRectangleShapeMinY = pointShape.Y;
                smallRectangleShapeMaxY = pointShape.Y + pixelHeight;
            }

            RectangleShape smallRectangleShape = new RectangleShape(smallRectangleShapeMinX, smallRectangleShapeMaxY, smallRectangleShapeMaxX, smallRectangleShapeMinY);

            GeoImage geoImage = rasterSource.GetImage(smallRectangleShape, 1, 1);
            Bitmap bitmap = new Bitmap(geoImage.GetImageStream(new GdiPlusGeoCanvas()));
            
            Color color = bitmap.GetPixel(0, 0);
            bitmap.Dispose();
            geoImage.Dispose();

            getPixelValueResult.PixelValue = new GeoColor(color.A, color.R, color.G, color.B);
            getPixelValueResult.IsValid = true;
            return getPixelValueResult;
        }
    }
}

using ThinkGeo.MapSuite.Core;

namespace ThinkGeo.MapSuite.Core
{
    class GetPixelValueResult
    {
        private bool isValid;
        private GeoColor pixelValue;

        public GetPixelValueResult()
            :this(false, new GeoColor())
        { }

        public GetPixelValueResult(bool isValid, GeoColor pixelValue)
        {
            this.isValid = isValid;
            this.pixelValue = pixelValue;
        }

        public bool IsValid
        {
            get { return isValid; }
            set { isValid = value; }
        }

        public GeoColor PixelValue
        {
            get { return pixelValue; }
            set { pixelValue = value; }
        }
    }
}



[David11]

Thankyou very much for this... sorry for the delay in replying, I'm about to get back to this problem next week.

Some initial impressions, I'd be interested to see the performance of this - not sure what sort of a hit the 1x1 bitmap creation will have as I'll be firing off tens of thousands of requests - nevertheless - you've given me the idea to make this work - I'd probably want to export a complete bitmap based on grid size that I want to query, where each pixel would compare exactly to my target grid - I would guess that there would be come distortion (averaging) of the pixel colours - but this would probably be a benefit ( a crude form of spatial averaging). I'll give both a go and compare... I'll let you know the result.

I've got another question for querying a vector layer which I will post as a new topic.

cheers

David McClymont

DHM Environmental Software Engineering.

[David]

David, 
  
   I’m glad I was able to spark some ideas.  If the pixels you are getting are around the same area then I suggest that maybe we create another function where you first pass in the bounds of the area where you will search, maybe include some limit and then we query the raster for that area.  Next we allow you to get pixel values based in the region which should be very fast since we are querying just one bitmap.  Another option is to add support for this at the lowest levels around where we call the third party libraries like MrSid and ECW.  In this way it is the quickest possible way to get the pixel.  of course if you know a region then its better to know that and it may even be faster. 
  
 David