Finding maximum pixel value

1642
6
04-17-2011 11:36 PM
PriyankaMehta
New Contributor
Hi all,

     I am trying to search for maximum pixel value in a raster within a defined search buffer and also at each degree from the source point.

So basically, if the source point is x,y . Starting from x,y, I want it to search maximum pixel value in 0°, 2°, 4° etc. till it reaches the buffer extent,  and the coordinates where maximum value is found.

Right now I am doing that by looping it and going to each pixel, checking its pixel value, comparing it to last value and then assigning the maximum. The code works, but at snail's pace.

I tried using IRasterstatistics and also tried exploring skipfactorX and y... but can seem to work it out.

My current code is in vba..
and it goes something like :


  For iRowCount = 0 To pRasterLayer.RowCount - 1
        For iColCount = 0 To pRasterLayer.ColumnCount - 1
            
            x = pRaster2.ToMapX(iColCount)
            y = pRaster2.ToMapY(iRowCount)

  Xle = Xle + 1
 
arrX(Xle) = x
arrY(Xle) = y

Next  iColCount
Next iRowCount

-------------------------------------------------------------
for q = 0 to xle

for dangle = 0 to 359 step 2

  
 ExtentX = 200 * Sin(rAngle)   ''200 is the search buffer distance. I am actually taking it from the user defined text box value
 ExtentY = 200 * Cos(rAngle)

    ExtentX = arrX(q) + ExtentX
    ExtentY = arrY(q) + ExtentY

  x2 = arrX(q)
 y2 = arrY(q)

 do until x2 = extentx and y2 = extenty

   ---- look through each pixel (Ipixelblock)
 ---- compare
---- assign max

loop

next
next




Another limitation with this code is that I have used arrays.
For a large raster and fine resolution, the size of the arrays exceeds to the limit where it gives me

runtime error : out of memory.
So I am stuck with using coarse resolution raster.




Please help

Thanks,
Priyanka
0 Kudos
6 Replies
BrianBottoms
New Contributor
If you have spatial analyst you might look at IZonalOp and ZonalStatisticsAsTable. If understand the �??degree from source point�?? part of the question, you would:

1. create a featureclass of all the lines emanating from the source point with a field containing the degree value.

2. Call something like �??

(dot net code�?�)
Dim zop As ESRI.ArcGIS.SpatialAnalyst.IZonalOp = _
New ESRI.ArcGIS.SpatialAnalyst.RasterZonalOpClass
Dim pRasterDataset As IRasterDataset
Dim pTable As ITable

pRasterDataset = pWs.OpenRasterDataset("MyRaster")
pTable = zop.ZonalStatisticsAsTable(CType(pFc, IGeoDataset), _
CType(pRasterDataset, IGeoDataset), True)

This will create an in memory table containing a max field. You can test the idea by just running ArcMap and digitizing a line over your raster then running the Zonal Statistics As Table tool and looking at the output. The difference when run from code is that you don�??t specify the zone field and instead the tool uses the featureclass object id to identify the line. So, you would need to loop the featureclass object ids, get the degree value then query the table for the matching ID and get the max value. I run this tool on very large grids with very large feature classes and it runs pretty fast.

Brian
0 Kudos
PriyankaMehta
New Contributor
hey thanks Brian,
  I just tried it manually. Haven't tried the code yet..
But i guess this should be it.

In my current code I am also retrieving the xy of the maximum pixel value.

do until x2 = extentx and y2 = extenty

   ---- look through each pixel (Ipixelblock)
 ---- compare
---- assign max
---- get xy of max pixel value

loop




Any ideas on how can I get the position (xy) of max pixel value by using zonal statistics.

Also, like I pointed before, I am using loops (for and do) for my purpose.

So basically, if I am in a situation where the line I have created has more than one pixel with the maximum value, the loop stops at the first maximum value and gives me the xy of that one. (It keeps iterating at that value till it reaches the extent)

Simplifying it, with my current code I am able to find out the first pixel of maximum value starting from the source (minimum distance) . Even if there are more pixels with the same value since I am looping it I can stop it at the first maximum value pixel and get its position. But then yes like I have already mentioned, it works extremely slowwwwww...


I hope I haven't confused you.

Can I have your email address . I can share the whole of my current code with you. That I guess would clear a lot of confusion on what I have and what I am looking for.

Mine is priyankamehta86@gmail.com




Thanks a lot..
Priyanka
0 Kudos
BrianBottoms
New Contributor
Ah, I misinterpreted your question. Statistics as Table will not work.

What you actually want is the closest point occurrence of the maximum value from a raster along a given azimuth (it is hard to phrase�?�).

Below is my approach. I use a structure to remember the x, y, azimuth, and raster value for each sample location from the start point to the search radius along a given azimuth. I then sort the collection of structures and pull the max value.

So, if I want to get all the locations for each even azimuth (2, 4, 6,�?�,360) at a search radius of 1000 meters with a sample point at every meter you call the function like:

Dim al as New ArrayList
al = MaxPoints(1000, 1, pPoint, pRaster)

This will return a collection of 180 points (pt structures), one for each even azimuth at the first occurance of the maximum value from the input raster. Runs pretty fast �?? hopefully faster than yours (if I understand the question now�?�).

Sorry if this just adds to the confusion. The code is a little sloppy. I guess if I wanted to be brief I'd say don't use pixel blocks and use pRaster.GetPixelValue(0, col, row)

Structure PT
        Implements IComparable
        Public X As Double
        Public Y As Double
        Public Degree As Integer
        Public Value As Double
        Public Position As Integer 'numeric order of point from center.

        Public Function CompareTo(ByVal obj As Object) As Integer Implements System.IComparable.CompareTo
            Dim tmpObj As PT = CType(obj, PT)
            Return (Me.Value.CompareTo(tmpObj.Value))
        End Function
    End Structure

Private Function MaxPoints(ByVal BufferRadius As Double, ByVal Interval As Integer, ByVal InPoint As IPoint, ByVal pRaster2 As IRaster2) As ArrayList

        'BufferRadius = search distance from center point (pPoint)
        'Interval = distance between sample points from center point to bufferradius
        '           should be more dense than raster cell size...
        'InPoint = center point
        'pRaster = raster with values for searching max value

        Dim al As New ArrayList ' temp array list to hold each sample point from center at given azimuth
        Dim finalAl As New ArrayList ' final return arraylist - always 180 points
        Dim col As Integer 'column of raster at sample point
        Dim row As Integer 'row of raster at sample point
        Dim newPT As PT ' structure to remember azimuth, sample point location, and raster cell value
        Dim x As Integer
        Dim j As Integer
        Dim z As Integer = 0
        Dim xCoord As Double
        Dim yCoord As Double
        Dim i As Double
        Dim count As Integer

        'Get max value for each even azimuth
        For x = 2 To 360 Step 2
            'Make sample point for every interval out to radius exenet along each even azimuth
            Do While z < BufferRadius
                xCoord = z * Math.Cos((90 - x) * Math.PI / 180.0F) + InPoint.X
                yCoord = z * Math.Sin((90 - x) * Math.PI / 180.0F) + InPoint.Y
                'Find corresponding raster cell at sample point location
                col = pRaster.ToPixelColumn(xCoord)
                row = pRaster.ToPixelRow(yCoord)
                'Get pixel value at sample point location.
                i = CDbl(pRaster.GetPixelValue(0, col, row))
                'Record sample point location, azimuth, and raster pixel value in PT structure.
                newPT.Degree = x
                newPT.X = xCoord
                newPT.Y = yCoord
                newPT.Value = i
                'Add to arraylist.
                al.Add(newPT)
                z += Interval
                count += 1
            Loop
            'Sort arraylist of PT structures on value Max value is now at the bottom
            al.Sort()
            'I know awkward...
            'Make sure that final location is the first occurrence i.e. get the last and compare to the next to last...
            newPT = CType(al(al.Count - 1), PT)
            For j = al.Count - 2 To 0 Step -1
                If CType(al(j), PT).Value = newPT.Value Then
                    If CType(al(j), PT).Position < newPT.Position Then
                        newPT = CType(al(j), PT)
                    End If
                Else : Exit For
                End If
            Next
            finalAl.Add(newPT)
            'Reset variables for next loop.
            al.Clear()
            z = 0
            count = 0
        Next
        'Final array list will have 180 PT structures - one for each even azimuth at the max cell location.
        Return finalAl
    End Function
0 Kudos
PriyankaMehta
New Contributor
Hey thanks again Brian,

   I ll try the code and respond asap.

Thanks a lot again !!
Regards,
Priyanka
0 Kudos
BrianBottoms
New Contributor
I cleaned up the code a little - I didn't know how to sort a collection of objects on 2 properties where the first is sorted descending (raster value) and the second sorted ascending (position of first max found). It turned out to be a lot easier than I thought. I was able to drop the more awkward part of the code. So the code does the same as before - returns 180 objects with x, y, and raster value of first maximum occurrence along each even azimuth from a source point (why you would do this I don't know...)

Structure PT
        Implements IComparable
        Public X As Double
        Public Y As Double
        Public Degree As Integer
        Public Value As Double
        Public Position As Integer 'numeric order of point from center.

        Public Function CompareTo(ByVal obj As Object) As Integer Implements System.IComparable.CompareTo
            Dim tmpObj As PT = CType(obj, PT)
            If Me.Value.CompareTo(tmpObj.Value) = 0 Then
                Return (Me.Position.CompareTo(tmpObj.Position))
            Else
                Return (tmpObj.Value.CompareTo(Me.Value))
            End If
        End Function
    End Structure
Private Function MaxPoints(ByVal BufferRadius As Double, ByVal Interval As Integer, ByVal InPoint As IPoint, ByVal pRaster2 As IRaster2) As ArrayList
        Dim al As New ArrayList
        Dim finalAl As New ArrayList
        Dim col As Integer
        Dim row As Integer
        Dim newPT As PT
        Dim z As Integer = 0
        Dim xCoord As Double
        Dim yCoord As Double
        Dim count As Integer

        For x As Integer = 2 To 360 Step 2
            Do While z < BufferRadius
                xCoord = z * Math.Cos((90 - x) * Math.PI / 180.0F) + InPoint.X
                yCoord = z * Math.Sin((90 - x) * Math.PI / 180.0F) + InPoint.Y
                col = pRaster.ToPixelColumn(xCoord)
                row = pRaster.ToPixelRow(yCoord)
                newPT.Degree = x
                newPT.X = xCoord
                newPT.Y = yCoord
                newPT.Value = CDbl(pRaster.GetPixelValue(0, col, row))
                al.Add(newPT)
                z += Interval
                count += 1
            Loop
            al.Sort()
            finalAl.Add(al(0))
            al.Clear()
            z = 0
            count = 0
        Next

        Return finalAl
    End Function



Brian
0 Kudos
PriyankaMehta
New Contributor
Thanks a lot again Brian,

The code looks neat to me. I am yet to test it though..
The reason I want to do it to determine the sky view factor. It is a value that determines how much sky is obstructed. The raster pixel values are nothing but the height of the buildings.

If i am standing in the midst of tall buildings the amount of sky i would see is the sky view factor.
The first tallest building that comes in my sight is the one that would be taken into consideration (Reason for the first maximum pixel) If there are consecutive pixels of the same elevation, then they might just be the breadth of that building

It is calculated by a formula and the output from above code is one of the input of the formula.

It is an interesting thing to study. Usually done for temperature modeling

Thanks again,
Priyanka
0 Kudos