How can I calculate the analysis grid area?

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