Este artículo presentará cómo extender la funcionalidad de la fórmula shoelace, para calcular el área de polígonos simples que no yacen en dos dimensiones, sino en un plano tridimensional. Digamos que los puntos \(P_{i} = x_{i},\ y_{i},\ z_{i};\ i = 1, 2, ... n\) representan los vértices de un polígono en el espacio. La siguiente visualización, creada con library(plotly), demuestra como luciría tal polígono.

Esta visualización interactiva demuestra que los puntos seleccionados yacen en un mismo plano; sin embargo, existen infinitas superficies curvas que también contienen esos mismos puntos. Para propósito de este artículo, nos concentraremos únicamente en superficies planas. De esta manera, resolveremos el problema de hallar el área plana encerrada por varios puntos tridimensionales.

Un plano en el espacio se define completamente a través de dos elementos: el vector normal (perpendicular) al plano en cualquiera de sus puntos, y uno de esos puntos. Sean \(\overrightarrow{n} = (a,\ b,\ c)\) el vector normal; \(P_0 = x_0,\ y_0,\ z_0\) un punto conocido; y \(P = x,\ y,\ z\) otro punto cualquiera del plano. Entonces el plano quedará definido a través de la ecuación (1)

\[\overrightarrow{n} \cdot \overrightarrow{P_0 P} = 0\]

\[(a,\ b,\ c) \cdot (x - x_0,\ y - y_0,\ z - z_0) = 0\]

\[ \begin{equation} a (x - x_0) + b (y - y_0) + c (z - z_0) = 0 \tag{1} \end{equation} \]

El razonamiento en esta ecuación es que \(\overrightarrow{n}\) es perpendicular a \(\overrightarrow{P_0 P}\); por lo tanto, el producto punto de ambos vectores debe ser igual a cero. La ecuación (1) suele ser reescrita de la siguiente forma:

\[a x + b y + c z = d;\ d = a x_0 + b y_0 + c z_0\]

Pero nosotros escribiremos la variable \(z\) como una función de \(x,\ y\). Más tarde se verá el por qué.

\[ \begin{equation} z = f( x,\ y ) = A x + B y + C;\ A = -\frac{a}{c};\ B = -\frac{b}{c};\ C = \frac{d}{c} \tag{2} \end{equation} \]

Regresando al planteamiento inicial, aunque no se conoce el vector normal \(\overrightarrow{n}\), es posible calcularlo empleando tres vértices del polígono. En efecto, para definir un plano tridimensional se necesitan al menos tres puntos (de la misma manera que, para definir una recta bidimensional, se necesitan al menos dos). Digamos que la siguiente matriz M contiene tres de los puntos representados en la visualización de arriba.

##     x  y z
## P1  8  2 6
## P2 -4 -2 7
## P3  5  0 5

Dado que todos se encuentran en el mismo plano, da lo mismo cuáles tres puntos son seleccionados. Las funciones points_normal() y points_plane() fueron escritas para calcular los componentes de \(\overrightarrow{n}\), y los coeficientes de la ecuación (2) respectivamente. Ambas funciones contienen operaciones algebraicas sencillas, y pueden ser revisadas en sus respectivos enlaces.

points_normal(M)
##    a    b    c
##  0.4 -1.0  0.8
points_plane(M)
## z = Ax + By + C
##     A     B     C
## -0.50  1.25  7.50

El resultado de points_plane(M) indica que el plano buscado es \(z = -\frac{1}{2} x +\frac{5}{4} y +\frac{15}{2}\). Ahora que conocemos su ecuación, podemos añadirlo a la visualización plotly.

También se ha añadido la proyección del polígono en el plano \(z = 0\), o su “sombra” por así decirlo. A esta sombra la denominaremos la región \(R\). Ahora, la razón de representar \(R\), y de haber escrito el plano como una función, es que se aplicará Cálculo Multivariado para resolver nuestro problema. De acuerdo al Cálculo, si existe una superficie determinada por \(z = f(x,\ y)\), el área de la porción de esa superficie correspondiente a los puntos de cierta región \(R\) viene dado por:

\[ \begin{equation} S = \iint\limits_R \sqrt{1 + \left( \frac{\partial f}{\partial x} \right)^2 + \left( \frac{\partial f}{\partial y} \right)^2} dS \tag{3} \end{equation} \]

\(S\) es el área superficial buscada y \(dS\) -el elemento diferencial de área- es igual a \(dxdy\) cuando las coordenadas son rectangulares. En el ejemplo ya planteado, \(R\) es este polígono bidimensional:

Normalmente es importante definir los límites de \(R\), a través de intervalos o funciones en \(x,\ y\). Además, es necesario calcular las derivadas parciales de \(f(x,\ y)\); si se tratase de una superficie curva, se obtendrá una expresión más o menos compleja que deberá ser integrada. Pero como se trata de un plano, con ecuación de la forma (2), las derivadas serán constantes: \(\frac{\partial z}{\partial x} = A;\) \(\frac{\partial z}{\partial y} = B\). Por ende no afectarán a la integral (3).

\[ \begin{equation} S = \sqrt{1 + A^2 + B^2} \iint\limits_R dS = \sqrt{1 + A^2 + B^2} \ S_R \tag{4} \end{equation} \]

Lo más importante es que nada será integrado, salvo el elemento diferencial de área. Esta integral equivale a calcular \(S_R\), el área plana de la región \(R\). Es un hallazgo útil pues, en vez de preocuparnos por los límites de \(R\), simplemente hay que calcular su área. Resumiendo: el área de un polígono que yace en un plano tridimensional es igual al área de su sombra (el mismo polígono pero con \(z = 0\)) multiplicada por cierta constante (dependiente de la ecuación del plano).

La fórmula (4) fue implementada en la siguiente función. Se trata de una extensión de la fórmula shoelace (la definición de shoelace_area() está en el artículo anterior) que multiplica su resultado por la constante adecuada, considerando un polígono en el espacio tridimensional.

shoelace_plane_area = function(vertices){
   suppressMessages(f <- points_plane(vertices))
   f = sqrt(1 + f["A"] ** 2 + f["B"] ** 2)
   names(f) = NULL; f * shoelace_area(vertices)
}

Resolviendo el ejemplo, shoelace_area(M) devuelve 22, mientras que shoelace_plane_area(M) devuelve 36.8951. La razón entre ambos valores es 1.6771, o también \(\sqrt{1 + \left( -\frac{1}{2} \right)^2 + \left( \frac{5}{4} \right)^2} = \frac{3}{4} \sqrt{5}\) si se expresa la razón con los coeficientes del plano.

Ahora plantearemos un ejemplo con datos más o menos reales. Se utilizará el dataset volcano que contiene la elevación -digitada a partir de un mapa topográfico- del volcán Maunga Whau.

Anteriormente se habló de cómo las funciones shoelace_area() y sf::st_area() permiten calcular lo mismo. Ahora se verá que también para shoelace_plane_area() ya existe una función equivalente: sp::surfaceArea(), que implementa un método (Jenness 2004) para calcular el área superficial -en lugar del área planimétrica- de un terreno, a partir de su modelo digital de elevación. Entonces, ejecutando sp::surfaceArea(volcano, 10, 10) -donde diez es el tamaño de las celdas en metros- se obtendrá el área del Maunga Whau, incluyendo el efecto de su topografía; este valor es de 55.9321 hectáreas.

Para replicar el funcionamiento de esta función hay que ser un poco laborioso. Sin entrar en detalles, el proceso es: se define una función que, dado un cuadrado con valores de elevación en cada esquina, divide el cuadrado en cuatro triángulos iguales y calcula sus áreas con shoelace_plane_area(); luego se itera esta función en cada cuadrado posible de la matriz volcano; finalmente, se insertan una fila y una columna extras con el menor área posible (100 m2) por celda. El último paso se efectúa para corregir el hecho de que la elevación debería estar al centro de cada celda, no en las esquinas.

matrix_area = function(m, d = 1){

   m = c(0, d, m[1],   0, 0, m[2],   d, 0, m[4],   d, d, m[3],   d/2, d/2, mean(m)) %>%
      matrix(5, , TRUE, list(c("P1", "P2", "P1", "P2", "P3"), c("x", "y", "z")))

   a = shoelace_plane_area(m[-c(3, 4), ])
   a = shoelace_plane_area(m[-c(1, 4), ]) + a
   a = shoelace_plane_area(m[-c(1, 2), ]) + a
   a + shoelace_plane_area(m[c(4,1,5), ])
}

volcano_area = matrix(0, nrow(volcano), ncol(volcano))

for (i in 1:nrow(volcano)) {
   for (j in 1:ncol(volcano)) {
      volcano_area[i, j] = if(i == 1 | j == 1) 100 else
         matrix_area(volcano[{i - 1}:i, {j - 1}:j], 10)
   }
}

sum(volcano_area)
## [1] 558755

De esta manera se estimó, con nuestras propias fórmulas, que la superficie del volcán es de 55.8755 hectáreas. Esto es el 99.8988 % del valor considerado más real. En el futuro, si se trata de superficies topográficas, se utilizará la función sp dedicada. Pero si alguna vez se tratase de polígonos en el espacio, tal vez la solución aquí descrita será la ideal. Para culminar, la visualización plotly del área real -por celda de 100 m2 planos- del 🌋

Referencias

Jenness, Jeff S. 2004. “Calculating Landscape Surface Area from Digital Elevation Models.” Wildlife Society Bulletin 32 (3): 829–39. https://doi.org/10.2193/0091-7648(2004)032[0829:CLSAFD]2.0.CO;2.