Subscribe to this thread
Home - General / All posts - Raster standard deviation
vincent

1,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 !

StanNWT
173 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.

vincent

1,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.

tjhb

9,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).

tjhb

9,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.

vincent

1,859 post(s)
#12-Mar-20 12:36

I can't say I'm really surprised. Don't give up !

tjhb

9,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.

tjhb

9,159 post(s)
#13-Mar-20 03:25

P.s. learning everything from Adam of course. Slowly but surely!

vincent

1,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 filter

Why ?

tjhb

9,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.

adamw


8,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.

vincent

1,859 post(s)
#15-Mar-20 13:05

A circular filter of radius 31 is what I use generally.

dyalsjas120 post(s)
#15-Mar-20 15:57

I also have a model that uses larger raster masks. The current functions (min, max, etc.) enable it, however the radius of the two masks are (circular) 15 and 151 pixels, respectively when using a 1m DEM. The smaller mask works well, the larger mask takes hours/days. Additional mask functions like standard deviation, are always welcome.

I suspect that as LiDAR data collections continue to commoditize high resolution elevation data, the ability to scale current analytic models (like hydrology models) to work with 1 meter resolution DEMs (instead of prevalent 5m, to 30m DEMs) will become more interesting.

Manifold already enables much of this functionality exceptionally well, but analytic models built to evaluate surface relationships using 5m to 30m DEMs lose a lot of performance/relevance when applied to a 1m DEM.

The first question from some of the other vendors I've talked with has been "Why not resample to a lower resolution?"

The second question tends to be "Can this be done in the cloud with multiple VM instances?"

The third question tends to be "Can you wait until we can build better multi-threading capability in "x" months/years/versions?"

Sure, all of these questions can be answered with "Yes", but Manifold usually disrupts that paradigm.

adamw


8,971 post(s)
#05-Apr-20 16:35

We've added TileFilterStDev in 9.0.171.2. There is a GPGPU variant, but due to memory constraints on the device, it will only work up to and including the radius of 8, calls with higher radius values will automatically be registered to a CPU. The CPU version, however, isn't bad performance-wise either, it uses single-pass robust math, etc. Try it out.

vincent

1,859 post(s)
#24-Mar-20 13:21

Any update Tim ?

tjhb

9,159 post(s)
#26-Mar-20 01:31

Like many others, we have just transitioned into full lockdown in New Zealand. (We are really really lucky here.) So now there is more time.

I admit to having been more than a bit shocked by your and dyalslas' ambitions for the radius, for me they caused a slowdown. And while I agree wholeheartedly with using a circular model (what else makes sense?), it does make things a little bit more tricky.

Maybe what I can contribute will be too basic to help very much.

Or, maybe, the idea of calculating SD over a radius of 151px--or even perhaps 31px--is, well, what exactly? In other words why. I suppose there could be a good reason, sometimes.

So what are the actual questions you are asking, in real terms (environmental terms, let's say). Why are you asking them?

It is really dumb (could even be idiotic) to throw the maximum currently available computing power at a poorly defined problem.

Analysis must be tailored to the question, question must closely follow the purpose.

Otherwise it is just using as many watts as current technology can consume. Which would be stupid. Define your problem.

But code is still forthcoming, yes. Just not on a dumb scale.

vincent

1,859 post(s)
#26-Mar-20 12:09

Hi Tim. Lockdown here too. But stil have a job and working from home. Things are well under control in health area.

It can be square, if it simplify the whole thing. No big deal.

All my project is individual-tree centered, from LiDAR data. Pixels are 1m width. 31m is the surrounding area of each tree affecting the tree (equal to max height of a tree). It would be larger from this standard deviation criteria, but 31 is the maximum allowed for CUDA in M8.

Take care !

vincent

1,859 post(s)
#26-Mar-20 22:51

It would be larger from for this standard deviation criteria, but 31 is the maximum allowed for CUDA in M8.

dyalsjas120 post(s)
#26-Mar-20 15:25

Tim,

I've re-read the documentation of the raster model I'm trying to recreate in Manifold. The model evaluates categories of surface variablility within a physical distance:

90 meter DEMs: 3 pixel radius sampling kernel

30 meter DEMs: 10 pixel radius sampling kernel

10 meter DEMs: 31 pixel radius sampling kernel

5 meter DEMs: 63 pixel radius sampling kernel

2 meter DEMs: 151 pixel radius sampling kernel

1 meter DEMS: 301 pixel radius sampling kernel.

For vertical relief variability, the maximum elevation within the search radius minus the minimum elevation within the search radius.

Slope diversity within the search radius

Aspect diversity (Aspect recoded into 17 categories, with 0 representing no signficant slope or aspect).

The key factor is that the model is based on a physical distances (approximately 300m).

Obviously, I could resample the data to a smaller pixel value (2, 3, or 5 meters), but I'd like to be able to evaluate the impacts of micro-relief (1 to 3 meters horizontal) that would be lost when resampling high-res lidar.

Manifold already has most of the functions built-in, but the data set I'm trying to evaluate is quite dense.

I started with (and still have a 1m DEM of Garrett County, Maryland that Art Lembo graciously shared with me in relation to this posting: http://www.georeference.org/forum/t142200.26#142201

It seemed ambitious to ask for 301 pixel radius masks, so I stepped back to 151 pixel radius, expecting to resample the DEM to 2m pixels.

I haven't tackled sending TileValues to a table and creating a layer of points, but given the underlying capability of Manifold to handle big data, I wonder if it would be more attainable to tackle the model that way, using buffers.

Mike Pelletier


1,698 post(s)
#26-Mar-20 17:41

Here is a old post regarding moving tile values to points to areas. It might be helpful for your project.

Dimitri


5,892 post(s)
online
#27-Mar-20 06:35

Another very easy way is using the transform, as shown in the Transfer DEM Heights to Areas in a Drawing example. Works for points too.

Manifold User Community Use Agreement Copyright (C) 2007-2019 Manifold Software Limited. All rights reserved.