# Finding maximum pixel value

968
6
04-17-2011 11:36 PM
Highlighted
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.

Thanks,
Priyanka
Tags (5)
6 Replies
Highlighted
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
Highlighted
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 valueloop`

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
Highlighted
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
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
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
'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
Highlighted
New Contributor
Hey thanks again Brian,

I ll try the code and respond asap.

Thanks a lot again !!
Regards,
Priyanka
Highlighted
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 StructurePrivate 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
Highlighted
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 