How would I go about calculating the area of the analysis grid that is above a given threshold value. I know I can find the percentage of all points above it, but I need the actual floor-area and preferably very accurately.

This Ecotect script demonstrates how to iterate over the current analysis grid in your model and calculate the actual surface area that is above a given threshold, in this case the current minimum scale value. Whilst Ecotect can already show the *percentage* of grid points above a threshold, this script includes the fractional areas of cells that are partially above and partially below.

A major part of this script is the interpolation of the effective surface area of those grid cells in which some of the node values are above the threshold and some below. Whilst it is possible to manually adjust the spacing of grid cells along each axis, such that individual grid cells can be of significantly different sizes, the analysis grid is always rectilinear so that the calculation of intersection areas is relatively straightforward. Figure 1 shows some examples of the equations required to determine the area of grid cells with different threshold intersections.

The right-most example shown in Figure 1 above shows that, in some instances, it is possible for the same threshold contour to pass through the same grid square twice. This effect is relatively easy to calculate, however it can occur on any combination of the four opposing corners, which makes it laborious to constantly test for.

Having always been a fan of the ''simpler is better'' principle, the code required to properly detect and handle all of these possible variations over the four edges of each cell was starting to balloon out pretty badly. Looking for a much simpler option that still demonstrated the same principles, I considered testing each cell as two separate triangles.

With a triangle, a particular contour line can only pass through it once. If it does, then there will always be either one corner above the threhold and two below, or two above and one below. This greatly simplifies the number of options that have to be checked for and handled.

Also, because we are only interested in the 2D area of the grid in order to derive the floor area (in this particular script), we only need to deal with the X and Y components of each grid cell's corner position. This makes the calculation of the effective surface area of a triangle pretty simple. Based on the triangle area equation given at the Math Open Reference website, we can start with a very simple function in our script.

-- Return the area of a triangle comprising 3x2D points. -- Dividing by 2x1000 converts millimetre positions into metres squared areas. function areaOfTriangle(Ax, Ay, Bx, By, Cx, Cy) return (abs((Ax * (By - Cy)) + (Bx * (Cy - Ay)) + (Cx * (Ay - By))) / 2000.0) end

We also need a couple of extra helper functions; the main one being for linearly interpolating the contour intersection value between two corners that cross the threshold.

-- A simple boolean test. function isAbove(value, threshold) if value >= threshold then return 1 end return 0 end -- Performs a linear interpolation between -- (v1 -> v2) given the ratio of u to (u1 -> u2). function Interpolate(u, u1, u2, v1, v2) -- If division by zero, return average of v1 and v2. if u2 ~= u1 then return((((u - u1) / (u2 - u1)) * (v2 - v1)) + v1) else return((v1 + v2) / 2.0) end end

Once we have these, we can create a function to test each triangle within a given grid cell and return its area. In Ecotect it is possible to hide parts of the analysis grid, so the first step is to check that all three triangle points are visible. Then. if all points are above the threshold, its entire area is returned. If all points are below, it returns zero area. We then need to check for either one or two points below the threshold.

-- ================================================ -- CELL ANALYSIS FUNCTIONS. -- rows(Y) -- di---C---ci 4---C---3 -- | / | | / | -- D / B D / B -- | / | | / | -- ai---A---bi cols(X) 1---A---2 -- ================================================ function getAreaOfTriangleAboveThreshold(threshold, index, cell) -- Check which trianglular facet to use. if index > 0 then ai = 4 else ai = 2 end totalArea = areaOfTriangle(cell.x[1], cell.y[1], cell.x[ai], cell.y[ai], cell.x[3], cell.y[3]) -- Check all nodes are visible. if cell.state[1] < 0 or cell.state[ai] < 0 or cell.state[3] < 0 then return 0.0 end -- Check if all nodes are above threshold. if cell.above[1] > 0.5 and cell.above[ai] > 0.5 and cell.above[3] > 0.5 then return totalArea end -- Check if all nodes are below threshold. if cell.above[1] < 0.5 and cell.above[ai] < 0.5 and cell.above[3] < 0.5 then return 0.0 end -- Cycle nodes. for ii = 1,3 do ai = ii -- Check to wrap next index. if ai >= 3 then bi = 1 else bi = ai+1 end -- Check to wrap next plus one index. if bi >= 3 then ci = 1 else ci = bi+1 end -- Check to invert triangle. if index > 0 then if ai == 2 then ai = 4 end if bi == 2 then bi = 4 end if ci == 2 then ci = 4 end end -- Check if previous and next segments rise above the threshold. if cell.above[ai] > 0.5 and cell.above[bi] < 0.5 and cell.above[ci] > 0.5 then Ax = Interpolate(threshold, cell.value[bi], cell.value[ai], cell.x[bi], cell.x[ai]) Ay = Interpolate(threshold, cell.value[bi], cell.value[ai], cell.y[bi], cell.y[ai]) Cx = Interpolate(threshold, cell.value[bi], cell.value[ci], cell.x[bi], cell.x[ci]) Cy = Interpolate(threshold, cell.value[bi], cell.value[ci], cell.y[bi], cell.y[ci]) area = areaOfTriangle(Ax, Ay, cell.x[bi], cell.y[bi], Cx, Cy) return totalArea - area end -- Check if previous and next segments dip below the threshold. if cell.above[ai] < 0.5 and cell.above[bi] > 0.5 and cell.above[ci] < 0.5 then Ax = Interpolate(threshold, cell.value[bi], cell.value[ai], cell.x[bi], cell.x[ai]) Ay = Interpolate(threshold, cell.value[bi], cell.value[ai], cell.y[bi], cell.y[ai]) Cx = Interpolate(threshold, cell.value[bi], cell.value[ci], cell.x[bi], cell.x[ci]) Cy = Interpolate(threshold, cell.value[bi], cell.value[ci], cell.y[bi], cell.y[ci]) area = areaOfTriangle(Ax, Ay, cell.x[bi], cell.y[bi], Cx, Cy) return area end end return 0.0 end

Finally we need to cycle through the grid and test each cell. As the overhead of continually requesting data from Ecotect is probably a bit greater than simply transferring some values, I chose to shuffle the data along as the loop cycles through the *j* and *k* directions in the grid.

-- ================================================ -- CYCLE THROUGH GRID. -- ================================================ areaAbove = 0 areaTotal = 0 for j = 1, grid.rows-1 do cell.value[1] = get("grid.cell", 0, j-1) cell.above[1] = isAbove(cell.value[1], grid.min) cell.x[1], cell.y[1] = get("grid.position", 0, j-1) cell.state[1] = get("grid.state", 0, j-1) cell.value[4] = get("grid.cell", 0, j) cell.above[4] = isAbove(cell.value[4], grid.min) cell.x[4], cell.y[4] = get("grid.position", 0, j) cell.state[4] = get("grid.state", 0, j) for i = 1, grid.cols-1 do cell.value[2] = get("grid.cell", i, j-1) cell.above[2] = isAbove(cell.value[2], grid.min) cell.x[2], cell.y[2] = get("grid.position", i, j-1) cell.state[2] = get("grid.state", i, j-1) cell.value[3] = get("grid.cell", i, j) cell.above[3] = isAbove(cell.value[3], grid.min) cell.x[3], cell.y[3] = get("grid.position", i, j) cell.state[3] = get("grid.state", i, j) -- Get cell square size. dx = cell.x[2] - cell.x[1] dy = cell.y[4] - cell.y[1] cellAbove = getAreaOfTriangleAboveThreshold(grid.min, 0, cell) cellAbove = cellAbove + getAreaOfTriangleAboveThreshold(grid.min, 1, cell) areaAbove = areaAbove + cellAbove cellTotal = (dx * dy) / 1000.0 areaTotal = areaTotal + cellTotal -- Shuffle values along. cell.value[1] = cell.value[2]; cell.above[1] = cell.above[2]; cell.state[1] = cell.state[2]; cell.x[1] = cell.x[2]; cell.y[1] = cell.y[2]; cell.value[4] = cell.value[3]; cell.above[4] = cell.above[3]; cell.state[4] = cell.state[3]; cell.x[4] = cell.x[3]; cell.y[4] = cell.y[3]; end end

The results are stored in the variables `areaAbove` and `areaTotal`, defined in lines 5 and 6 above. We can then print these out or display them as part of a more complex report.

Attachment | Size |
---|---|

GetGridArea.zip | 2.17 KB |