Home - General / All posts - Raster standard deviation
 vincent1,859 post(s) #11-Mar-20 13:33 Hi, I'm looking for a way to obtain the standard deviation for each pixel of an image based on its neighborhood. Preferably in M8, but M9 is fine too. My current method calculates standard deviation for a vector grid cell of 1 hectare with : Update(Select * from(SELECT [a].[ID], STDEV([Height (I)]) as [res]FROM [Fishnet_1ha_inters_planif]as [a],[MHC_pg10]WHERE [MHC_pg10].[Invisible (I)] = False and Touches([a].[ID],AssignCoordSys(NewPoint([MHC_pg10].[Center X (I)],[MHC_pg10].[Center Y (I)]) ,COORDSYS ("MHC_pg10" as COMPONENT)))GROUP BY [a].[ID] ) as [D2] join [Fishnet_1ha_inters_planif] as [b] on [b].[ID] = [D2].[ID])Set [std_dev] = [res]I would rather like to obtain the results based on a moving window for each cell. I didn't found any function to do that either in 8 or 9.Thank you !
 StanNWT173 post(s) #11-Mar-20 15:21 This is essentially a matrix filter of a user defined size where the calculation is the standard deviation correct? I bet ArcGIS Desktop or ArcGIS Pro would crash, given it's inability to do anything other than a 3 x 3 matrix filter and certainly not a standard deviation on that matrix. Not that they couldn't try to program it in python.
 vincent1,859 post(s) #11-Mar-20 15:54 This is essentially a matrix filter of a user defined size where the calculation is the standard deviation correct?Yes.
 tjhb9,159 post(s) #11-Mar-20 15:59 It's pretty straightforward in 9 or 8. (It will be interesting to see the speed difference--how much CUDA processing has been improved in 9.)I'll post code after breakfast (or coffee).
 tjhb9,159 post(s) #12-Mar-20 07:01 After I woke up properly I found it was less straightforward, no surprise to you I suppose.More or less sorted now, code tomorrow.
 vincent1,859 post(s) #12-Mar-20 12:36 I can't say I'm really surprised. Don't give up !
 tjhb9,159 post(s) #13-Mar-20 03:09 Won't give up! Bear with me a bit longer if you can.I'm working on two different approaches, both great learning experiences (and good examples).A summary for your patience...One approach uses TileCutRect() to define the analysis window around every pixel. (So the new tiles systematically overlap, by the window radius. On the other hand each tile is small.) These window tiles are passed to a custom aggregate function to perform the standard deviation.The second approach uses a stack of custom matrix filters. (We need a separate filter for each neighbouring pixel, since their differences from the average need to be independently squared. Squaring can't be implented within a matrix filter, although it is an efficient tile operation in itself.)I expect the second approach to be significantly faster, since I think it can be performed entirely on GPU. I don't think that can be done for the first approach.More soon.
 tjhb9,159 post(s) #13-Mar-20 03:25 P.s. learning everything from Adam of course. Slowly but surely!
 vincent1,859 post(s) #13-Mar-20 18:13 Bear with me a bit longer if you can.I can wait 2 weeks if needed.Squaring can't be implented within a matrix filterWhy ?
 tjhb9,159 post(s) #13-Mar-20 19:02 Custom matrix filters allow us to define arbitrary coefficients for a window including a centre pixel and its neighbours (using the functions with Def in their name), then to apply a mathematical operation over the source pixels in each window, using the coefficients as weights.The operations we have currently (the functions without Def) are average (plain TileFilter, also used for blurring, smoothing, sharpening, edge detection), sum, min, max, median, major, diversity. These all inherently involve neighbours.Not multiplication, division, powers or roots. Those operations are available (with many others) for each individual pixel in a tile (while processing a tile at once), but not for matrix filters that operate amongst neighbours.In any case the squaring needed here is on individual pixels (the results of separate matrix filters), not a relation between each pixel and one or more neighbours.
 adamw8,971 post(s) #15-Mar-20 11:20 Missed this thread somehow.It struck me that while you can solve the problem with TileCutRect - which is what I would do, because that's simplest, although not terribly fast - it would have been much better if in addition to TileFilterMax / Min / Sum, etc, we simply had TileFilterStDev (+ variants). I posted a request to add that into the wishlist. :-)vincent, how big of a window you are planning to use around each pixel? 5x5, 7x7, something like that? Because if it is more like 101x101, be prepared to either reduce the resolution first (to slash the number of pixels and the size of the window - because I presume you want to take all pixels within X meters and reducing the resolution will reduce the number of pixels in that X) or be waiting for a fair amount of time, regardless of the specific solution. Processing a 50 MB raster with a 101x101 window applied to each pixel means having to process 50 MB * 10,000 = 500 GB of data, this is going to take some time no matter what.
 vincent1,859 post(s) #15-Mar-20 13:05 A circular filter of radius 31 is what I use generally.