# How can I calculate the analysis grid area?

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.

Figure 1 - Some examples of the various contour intersections with a rectangular grid cell.

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.

Figure 2 - The four variations of contour intersections with a triangular facet.

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
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)

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.

AttachmentSize
GetGridArea.zip2.17 KB

### Authors/Reviewers

View desktop or mobile version of site.