Un índice de rugosidad del terreno a escala municipal –updated–

1. Introducción

Este trabajo es una mejora y actualización de Goerlich y Cantarino (2010a). En aquella ocasión, hace ya más de 10 años, elaboramos un índice de rugosidad del terreno a escala municipal¹ a partir de Modelos Digitales de Elevaciones (MDE) públicos y con cobertura mundial. ¡El Centro de Descargas del Centro Nacional de Información Geográfica (CNIG), dependiente del Instituto Geográfico Nacional (IGN), no existía por aquel entonces!

¡Mucho ha ☔ desde entonces! y disponemos en la actualidad de información geográfica abundante, libre y de excelente calidad que permite mejorar y actualizar lo que hicimos en su momento.

En particular, Goerlich (2022), a partir del Modelo Digital del Terreno (MDT)² con paso de malla de 5m –celdas de 25m²– disponible en el Centro de Descargas del CNIG, elabora un modelo digital, con la misma resolución, en el que el valor de cada celda representa la superficie del paisaje de la misma. Es decir, si en dicho modelo las celdas cubren una superficie plana, 2D, sobre el terreno de 25m², el valor del pixel en dicho modelo constituye una estimación de la superficie teniendo en cuenta la orografía del paisaje, es decir, la rugosidad del terreno. Lo que denominaremos la superficie 3D.³

La ratio de la superficie 3D y la superficie 2D de una determinada área constituye un índice sencillo e interpretable de rugosidad, carente de unidades de medida y tremendamente intuitivo. Aunque, ciertamente, índices de rugosidad hay muchos, y normalmente no requieren una estimación de la superficie 3D, siempre parten de las pendientes, y el MDT con paso de malla de 5m del IGN constituye una información excelente como punto de partida.

Así pues, se trata de una mejora y actualización de Goerlich y Cantarino (2010a) en todas las direcciones. Cambiamos la información de base, un modelo nacional de alta resolución frente a modelos globales con menor resolución, y cambiamos de índice, frente a uno de los más populares (Riley, DeGloria y Elliot 1999), con unidades de medida, y por tanto dependiente de la resolución del modelo de partida, a uno más intuitivo e interpretable, y que además tiene una agregación consistente, 👍.

La estructura del trabajo es la siguiente. El apartado 2 describe el método seguido para la estimación de la superficie 3D de cada una de las celdas del MDT y el siguiente relaciona la información utilizada. El apartado 4 examina muy brevemente cuestiones metodológicas relacionadas con la superficie 3D y la rugosidad del terreno, así como las propiedades de los índices utilizados. El apartado 5 muestra los resultados obtenidos a nivel de municipio, contraponiendo altitud frente a rugosidad, y el 6 relaciona los resultados anteriores con la distribución de la población en el largo plazo. Un último apartado concluye.

2. Estimación de la superficie del paisaje o 3D

Si la medición de la superficie 2D de un área extensa, como una región o un país, no parece tarea sencilla, la estimación del paisaje –3D– se antoja un trabajo mucho más arduo (Goerlich 2022). Partiendo de un MDT de elevada resolución realizaremos estimaciones de la superficie 2D y 3D de los municipios, para construir un índice de rugosidad a partir de ellas.

Aunque es cierto que podemos encontrar las superficies municipales 2D en el Centro de Descargas del CNIG –Nomenclátor Geográfico de Municipios y Entidades de Población (NGMEP)– preferimos realizar nuestras propias estimaciones a partir de la misma fuente de información que utilizamos para construir las estimaciones de la superficie 3D. Ello proporciona una mayor homogeneidad y consistencia en los resultados, al fin y al cabo nuestro interés aquí se centra en la ratio entre superficies y parece lógico que ambas utilicen la misma información de partida y los mismos métodos.

El MDT utilizado es, en realidad, un modelo raster en coordenadas proyectadas. Dada su resolución –25m² por celda–, y el hecho de que todas las celdas tienen el mismo tamaño la superficie 2D consiste en contar el número de pixeles dentro del polígono de referencia, y multiplicar dicho resultado por 25m². Utilizando este procedimiento, y la información descrita en el apartado siguiente, Goerlich (2022) ofrece una estimación de la superficie 2D de España de 506,239.59 km².

El método utilizado para la estimación de la superficie 3D se debe a Jennes (2004) y consiste, básicamente, en construir triángulos 3D partiendo de las elevaciones de un MDT, a partir de los cuales podemos determinar fácilmente la superficie 3D.

Con la información de un MDT –celdas cuadradas y todas ellas del mismo tamaño– en el que el valor del pixel represente la altitud del terreno procedemos de la siguiente forma. Tomando una celda central de referencia y sus 8 vecinas –queen contiguity–, los valores de las elevaciones de todas estas celdas –9 en total– pueden utilizarse para calcular las diferencias de altura entre cada par de celdas vecinas. Estas diferencias, junto con la resolución del modelo raster, el Teorema de Pitágoras y trigonometría elemental –que simplemente utiliza el coseno del ángulo de las pendientes–, permiten determinar la longitud de todos los bordes de los 8 triángulos levantados en 3D formados por la conexión entre el centro de la celda de referencia y el centro de las 8 celdas vecinas. Dadas las longitudes de todos estos bordes –16 en total–, convenientemente recortados por la base de la celda de referencia, podemos estimar el área de cada uno de los 8 triángulos por medio de la fórmula de Heron (Abramowitz y Stegun 1972). La suma del área de todos estos triángulos –8 en total– es una estimación de la superficie 3D de la celda de referencia. Los dos paneles de la Figura 1 ilustran el procedimiento.

---

---

El método se basa, pues, en deformar la superficie plana de la celda de referencia de acuerdo con las pendientes entre dicha celda y sus 8 vecinas. Esta deformación está constituida por 8 triángulos que se levantan en 3D y cuya área, sumada para todos estos triángulos, supone una estimación en 3D de la superficie de la celda de referencia. Pasamos de esta forma de un MDT que representa altitudes a un MTD cuyos valores representan la superficie del paisaje de la celda en cuestión. La resolución de la información obtenida es idéntica a la resolución de la información de partida.

La suma de los valores de este nuevo raster, recortado por el polígono de referencia, constituye una estimación de la superficie 3D para dicho polígono. Utilizando este procedimiento, y la información descrita en el apartado siguiente, Goerlich (2022) ofrece una estimación de la superficie 3D de España de 526,781.20 km². Ello representa un crecimiento, para el conjunto nacional, del 4.06%, lo que implica una pendiente media del 28.8% o de 16.1º sexagesimales.

Goerlich (2022) muestra claramente como los resultados son muy sensibles a la resolución de la información de partida y también, aunque en menor medida, al método de estimación. Este trabajo explota, a nivel municipal, la mayor resolución de la información elaborada en aquel trabajo con el método más preciso de los ensayados.

3. Descripción de la información utilizada

La información utilizada procede básicamente de dos fuentes. Un MDT con información sobre la altitud del terreno con elevada resolución y sobre el que implementar el procedimiento que acabamos de describir. Y una capa de contornos administrativos que delimite los objetos espaciales en los que estamos interesados, municipios, y por agregación acote el espacio para el que estimamos la superficie. Esta capa se utiliza para recortar –mask– la información raster.

4. Rugosidad del terreno y superficie 3D

La rugosidad del terreno puede ser definida de forma genérica como una característica relacionada con la irregularidad o accidentalidad del territorio. Sin embargo, una cuantificación de esta propiedad es, por su propia naturaleza, de difícil implementación. Para biólogos y naturalistas la medición de la rugosidad del terreno se ha convertido en una variable clave que ayuda a explicar los hábitats utilizados por especies de animales terrestres (Gionfriddo y Krausman 1986; Wiggers y Beasom 1986) o aspectos relacionados con la diversidad vegetal (Crawford 2008). No obstante esta característica del terreno también ha recibido cierta atención en geología (Hobson 1972) y en determinadas aplicaciones de ingeniería (Polidori, Chorowicz y Guillande 1991; Aguilar, Aguilar, Agüera, Carvajal y Sánchez 2002).

Este apartado ofrece un muy breve repaso de la literatura, explora como algunos autores han hecho operativa la idea genérica de rugosidad y muestra la relación entre la superficie 3D y la rugosidad.

En el campo del estudio de los hábitats de especies animales, el primer indicador de rugosidad del terreno se debe a Beasom, Wiggers y Giordono (1983), y se basa en el supuesto de que la rugosidad es función de la longitud total de las curvas de nivel de una determinada área. Propuesto antes de la explosión de la cartografía digital es un índice intuitivo, si bien su elaboración es intensiva en trabajo o solamente factible para determinadas áreas de tamaño reducido si su elaboración se realiza a partir de cartografía tradicional en papel. Quizá por lo temprano de la aportación, en relación a la información y tecnología disponible para su cálculo, lo cierto es que la propuesta de estos autores no ha alcanzado difusión posterior.

La aparición de los modelos cartográficos digitales favoreció la propuesta de otros índices de rugosidad que fueran sencillos de calcular mediante técnicas de Sistemas de Información Geográfica (GIS) y comparables entre áreas distantes y de diferente naturaleza. La primera de estas propuestas se debe a Riley, DeGloria y Elliot (1999), la idea es relativamente simple y es la propuesta que ha triunfado en las aplicaciones prácticas, aunque dispone de múltiples variantes.

A partir de un MDT, la metodología para obtener los valores de rugosidad a nivel de celda se basa en el análisis de las celdas vecinas, según se representa en el panel superior izquierdo de la figura 1. Para cada celda disponemos de un valor de altitud (referido al centro de la cuadrícula), de forma que es posible calcular un indicador de la rugosidad combinando las pendientes entre la celda central y todas sus adyacentes –sólo las celdas contiguas, 8 en total–.¹⁵ Un Índice de Rugosidad del Terreno (TRI) debe combinar de alguna forma todas las pendientes de las celdas adyacentes (8 en nuestro ejemplo). Es posible tomar como índice simplemente el valor máximo o mínimo de todas estas pendientes contiguas, pero un índice global debe situarse razonablemente entre estos dos valores extremos.¹⁶

Una posibilidad es simplemente promediar de algún modo estas pendientes o diferencias de altitudes (Wilson, O’Connell, Brown, Guinan y Grehan 2007), otra utilizar métodos algo más sofisticados, como la propuesta original de Riley, DeGloria y Elliot (1999), que consiste en tomar la raíz cuadrada de la suma de los cuadrados de las diferencias de altitud de cada celda vecina respecto a la celda central. Se trata de un índice muy parecido a una versión no escalada de la varianza, donde las desviaciones se miden respecto a la celda central. Se obtiene de esta forma un índice de rugosidad a nivel de celda, que deberá ser promediado de algún modo si lo que queremos es la rugosidad de áreas más extensas.

El índice de rugosidad de Riley, DeGloria y Elliot (1999), o alguna de sus variantes, ha hecho fortuna en el análisis empírico de esta característica del terreno, no solo porque es sencillo y poco costoso de implementar, sino porque la práctica totalidad del software dispone actualmente de comandos o scripts para su obtención, en muchos casos con la posibilidad de elegir entre varias alternativas.¹⁷ Este es el índice utilizado en Goerlich y Cantarino (2010a), quienes ofrecen intuición adicional sobre el mismo y examinan algunas de sus extensiones posibles, incluida la posibilidad de extenderlos en la dirección de la literatura de la desigualdad, ya que el índice de rugosidad de Riley, DeGloria y Elliot (1999), y sus variantes, es una medida de dispersión de las diferencias de altitud –elevaciones– entre las celdas vecinas respecto a la central. Dicho índice depende de las unidades de medida y de la resolución del raster de partida, y en consecuencia los intervalos de clasificación del terreno que se ofrecen en la literatura no son válidos con generalidad.

Es posible encontrar en la literatura otras propuestas sobre índices de rugosidad que no sean variaciones sencillas de la propuesta de Riley, DeGloria y Elliot (1999). Goerlich y Cantarino (2010a) mencionan dos.

En primer lugar, la de Sappington, Longshore y Thompson (2007), que adaptan ideas previas de Hobson (1972). Estos autores observan que el índice propuesto por Riley, DeGloria y Elliot (1999) está fuertemente correlacionado con la pendiente, y en este sentido es incapaz de distinguir entre terrenos con fuertes pendientes, pero relativamente poco rugosos, de terrenos que teniendo pendientes importantes sean irregulares y quebrados, es decir, tengan una rugosidad elevada. El argumento principal es que una medida de rugosidad debe incorporar la variabilidad en dos atributos que siempre es posible medir a partir de un MDT: (i) la pendiente (slope), y (ii) el aspecto (aspect).¹⁸ Atributos que pueden calcularse rutinariamente con los programas de GIS. Estos dos aspectos definen un plano para cada celda, de forma que la idea de Sappington, Longshore y Thompson (2007) es descomponer los vectores ortonormales a este plano para cada celda en sus tres componentes mediante el uso de trigonometría estándar, sumar dichas componentes por separado y calcular la norma euclídea del vector resultante. A partir de aquí proponen un índice de rugosidad vectorial que carece de unidades, y varía entre 0 –terreno totalmente plano o de mínima rugosidad– y 1 –máxima rugosidad–.

Aunque este índice tiene su atractivo,¹⁹ al ser la variabilidad en el aspecto de las celdas la que trata de recoger lo quebradizo del terreno, también presenta ciertos inconvenientes que no son resaltados por Sappington, Longshore y Thompson (2007). Goerlich y Cantarino (2010a) muestran como el tratamiento simétrico de la pendiente y el aspecto –dirección– que hace el índice puede conducir a valores de difícil interpretación. Por ejemplo, casos –extremos– en los que una celda tenga la máxima pendiente –90º–, y todas la celdas implicadas en el cálculo presentaran un aspecto también de 90º, produciría un valor del índice igual a 0, ¡es difícil pensar que estamos ante una situación de rugosidad nula!, ¡desde luego no es terreno plano!, ¡👎! Por otra parte, el límite superior a 1 del índice es un límite teórico, en la práctica son de esperar valores mucho más bajos.²⁰

Ciertamente Sappington, Longshore y Thompson (2007) están en lo correcto cuando argumentan que los índices de rugosidad convencionales, como el de Beasom, Wiggers y Giordono (1983) o el de Riley, DeGloria y Elliot (1999) y sus variantes, están muy correlacionados con la pendiente del terreno –de hecho están construidos a partir de esta característica– y que no recogen un aspecto relacionado con la rugosidad, pero independiente de la pendiente, como es lo quebradizo del paisaje. También parece cierto que esta es una cuestión relacionada con la resolución del MDT de partida, ya que un tamaño más pequeño de la celda debe ser capaz de reflejar mejor la rugosidad del terreno a gran escala.

En segundo lugar, la de Jennes (2004). La preocupación fundamental de este autor no está relacionada directamente con el desarrollo de una medida cuantitativa de la rugosidad del terreno, sino con la medición de la superficie de una determinada área teniendo en cuenta la orografía del terreno, es decir, la superficie 3D. De esta forma, a partir de un MDT, desarrolla el método de estimación de dicha superficie expuesto en el segundo apartado del trabajo.

Una vez disponemos de estimaciones de la superficie 2D y 3D –obtenidas con el mismo modelo de datos y la misma resolución–, transformar ambas estimaciones en un índice de rugosidad es trivial. La idea es relativamente sencilla. Si la superficie 2D de un área, sin tener en cuenta la orografía, son 100 km², y la superficie 3D son 110 km² la ratio de ambas superficies,

$$TRI_1 = \frac{Sup_{3D}}{Sup_{2D}}$$

que no es más que el crecimiento de la superficie al tomar en cuenta la orografía del terreno, constituye una medida de rugosidad del mismo. En nuestro ejemplo, la ratio es 1.1, lo que es equivalente a un crecimiento del 10% al pasar de la superficie plana a la superficie que toma en cuenta la orografía. Ello representa una pendiente media del 46% –notablemente elevada– o, de forma equivalente, un ángulo de 25º sexagesimales.

Como índice de rugosidad, la ratio $Sup_{3D}/Sup_{2D}$ carece de unidades de medida, su límite inferior es 1 –terreno plano– y no está acotado superiormente. Presenta la ventaja de tener una interpretación sencilla –el crecimiento en la superficie al tomar en cuenta la orografía del terreno– y ser muy intuitivo.

Si se quiere un índice acotado en el intervalo $[0, 1]$ podemos considerar la inversa del índice anterior, $Sup_{2D}/Sup_{3D}$, pero ahora un valor mayor del índice representa una menor rugosidad, por lo que si queremos la interpretación habitual –a mayor valor mayor rugosidad– resulta natural restar dicha ratio de la unidad,

$$TRI_2 = 1 - \frac{Sup_{2D}}{Sup_{3D}}$$

que tiene la interpretación de la pérdida relativa de superficie al pasar de la superficie 3D a la superficie 2D. Dicho índice está acotado en el intervalo $[0, 1]$, toma el valor 0 para terreno plano y se incrementa conforme aumenta la rugosidad del terreno. Aunque su valor máximo teórico es 1 este nunca se alcanzará, puesto que la superficie 3D nunca será infinita.²¹

Obsérvese que $TRI_1 - 1$ y $TRI_2$ simplemente escalan la misma variación absoluta de superficies, $Sup_{3D} - Sup_{2D}$, por una referencia diferente, $Sup_{2D}$ en el primer caso y $Sup_{3D}$ en el segundo.

Naturalmente ambos índices, $TRI_1$ y $TRI_2$, son ordinalmente equivalentes, aunque presentan una cardinalidad diferente. Desde el punto de vista de la cardinalidad $TRI_1$ es claramente preferido sobre $TRI_2$, ya que tiene una interpretación inequívoca.

En la práctica, todos los índices de rugosidad mencionados se calculan a nivel de celda, incluido el de Jennes (2004), ya que la superficie 3D se calcula para cada celda –tal y como se expone en el segundo apartado–. Para áreas más grandes, como municipios o provincias es necesario agregarlos por algún procedimiento, y aquí los índices basados en el cálculo de la superficie 3D presentan una clara ventaja, ya que permiten una agregación consistente.

Tanto los índices basados en las pendientes (Riley, DeGloria y Elliot 1999) como el propuesto por Sappington, Longshore y Thompson (2007) suelen obtenerse para áreas grandes mediante promedio de los valores de las celdas contenidas en dichas áreas.²² Los índices basados en el cálculo de superficies 2D y 3D permiten la agregación natural de sus componentes, y a partir de las superficies 2D y 3D agregadas se obtiene la ratio entre ambas. No es necesario –ni conveniente– agregar los índices de rugosidad a nivel de celda, sino proceder a la agregación de superficies y calcular el índice a partir de las superficies agregadas. Una propiedad de consistencia en la agregación que no es satisfecha por el resto de índices. Naturalmente esta ha sido la agregación utilizada en este trabajo a la hora de elaborar los índices de rugosidad a nivel municipal.²³

5. Resultados para los municipios españoles: Altitud versus rugosidad.

Este apartado compara la altitud con la rugosidad a nivel municipal, donde la rugosidad se obtiene como la ratio de las superficies 2D y 3D, y esta última ha sido calculada por el método de Jennes (2004) a partir del MDT05. En particular, utilizaremos el índice $TRI_1$, que es simplemente la ratio de superficies, $Sup_{3D}/Sup_{2D}$, por ser más intuitivo y tener una traducción directa a la pendiente media del área considerada.²⁴

Que la rugosidad constituye una característica del paisaje diferenciada de la altitud es evidente a partir de la inspección visual de los Mapas 1 y 2. El Mapa 1 muestra la altitud media a nivel municipal²⁵ y permite identificar claramente las cadenas montañosas, la meseta y las depresiones de nuestro país. ¡Es una imagen que nos es familiar de cuando estudiábamos geografía en el bachiller, 👍! Observamos una distribución relativamente suave de esta variable y permite identificar a España como un país de contrastes, montañoso en algunos lugares pero con amplias depresiones –Ebro y Guadalquivir– en otros.

---

---

El Mapa 2 muestra el índice de rugosidad $TRI_1$ calculado a partir del MDT05 y agregado de forma consistente a nivel municipal.²⁶ ¡Esta imagen nos es menos familiar, 👎! La distribución espacial de la rugosidad es mucho menos suave que la de la altitud, ¡al fin y al cabo la rugosidad viene a ser, aproximadamente, la primera derivada de la altitud y en consecuencia tiene un comportamiento más abrupto! Podemos observar amplias zonas costeras, de baja altitud, pero elevada rugosidad, y al mismo tiempo áreas elevadas que se corresponden con valores bajos de la rugosidad. Todo indica que ambas dimensiones del paisaje contienen información diferente sobre la orografía del terreno. La correlación entre ambas variables a nivel municipal es reducida, aunque ligeramente positiva, 0.177.²⁷

---

---

Las diferencias entre altitud y rugosidad a nivel municipal son evidentes si comparamos su distribución a esta escala –Figura 2–. Ambas distribuciones son asimétricas, con una larga cola superior, pero ello es mucho más evidente en la distribución de la rugosidad. La moda para la altitud se sitúa algo por encima tanto del valor medio como de la mediana, a unos 800m, mientras que la de la rugosidad se localiza muy por debajo de la media, incluso por debajo del primer cuartil, y relativamente cerca del valor mínimo del índice, 1. Dicha distribución está además mucho mas concentrada –obsérvese el eje de ordenadas en ambos histogramas que muestran el número de municipios y tienen escalas diferentes–, su coeficiente de variación es de 0.042, frente a 0.544 del coeficiente de variación de la distribución de la altitud. Lo mismo sucede con el índice de Gini, 0.020 para la rugosidad frente a 0.307 para la altitud.

---

---

Efectuando cuatro cortes de altitud, (i) hasta los 300m, (ii) más de 300 y hasta los 700m, (iii) más de 700 y hasta los 1000m y (iv) por encima de esta, de forma que tenemos aproximadamente dos cortes por debajo de la altitud media de España y dos por encima, obtenemos los resultados que aparecen en la Tabla 1.²⁸

Con esta partición, en cada extremo tenemos casi la quinta parte de los municipios –y un reparto similar de la superficie–, pero la distribución de la población es radicalmente diferente. Es conocido que la población se localiza mayoritariamente en la costa y en el valle (Goerlich y Molla 2021) y en consecuencia en el primer intervalo tenemos más de la mitad de la población, el 53.5%, pero en el último tan solo un 2.6%, 1.2 millones de personas de los más de 47 millones de residentes en España. El hecho más destacable de la Tabla 1 es, probablemente, la ausencia de correlación entre altitud y rugosidad si eliminamos los municipios por encima de los 1000m de altitud media, lo que podríamos identificar como municipios de montaña. No parece existir tendencia apreciable entre $TRI_1$ y altitud hasta que no sobrepasamos los 1000m, y esta es ya una altitud considerable. De hecho la correlación entre ambas variables excluyendo el último intervalo es prácticamente nula, –0.026–.

Tabla 1. Indicadores municipales por cortes altimétricos

	Municipios1		Población		Superficie		Rugosidad		Pendiente
Altitud	n	%	N	%	km²	%	Indice	Base 100	%	°
Hasta 300m	1,544	18.8	25,341,498	53.5	86,749	17.1	1.028	98.75	23.6	13.3
(300, 700]	2,377	28.9	16,060,316	33.9	175,787	34.7	1.037	99.67	27.5	15.4
(700, 1000]	2,804	34.1	4,748,079	10.0	150,577	29.7	1.034	99.38	26.4	14.8
Más de 1000m	1,487	18.1	1,235,214	2.6	93,126	18.4	1.070	102.78	37.9	20.8
España	8,212	100.0	47,385,107	100.0	506,240	100.0	1.041	100.00	28.8	16.1
1 El número de municipios incluye los territorios mancomunados
Fuente: MDT05 del IGN, Padrón 2021 del INE y elaboración propia

Ya hemos observado que la distribución de la rugosidad es tremendamente asimétrica a escala municipal –Figura 2–. El valor de $TRI_1$ para España es de 1.041, lo que equivale a una pendiente media del 28.8% –o de forma equivalente de 16.1º–. Es una pendiente media elevada, lo que ciertamente justifica la afirmación de que España es un país montañoso más allá de su altitud. No solo tiene una altitud media considerable dada su longitud de costa, sino que además es ciertamente abrupto en una gran parte de su geografía. El 66% de los municipios tienen un valor de $TRI_1$ inferior al promedio nacional, un 91% de los municipios lo tiene inferior al 1.1, un valor que representa un crecimiento del 10% al pasar de la superficie 2D a la 3D, o una pendiente media del 45.8%, y un 99% inferior al 1.2, lo que representa una pendiente media del 66.3%.

La Tabla 2 identifica los 25 municipios con mayor rugosidad.²⁹ Todos ellos con valores de $TRI_1$ por encima de 1.255, lo que supone una pendiente media superior al 75%, si bien el valor máximo, 1.392 alcanzado por San Andrés y Sauces –38033–, Santa Cruz de Tenerife, no alcanza el 100% de pendiente, lo que representaría un ángulo promedio en su término municipal de 45º, ¡aunque se queda cerca, 96.8%, –44.1º–! Curiosamente menos de la mitad de estos municipios –11– superan los 1000m de altitud promedio y la mayoría son pequeños, aunque encontramos una capital de provincia, Santa Cruz de Tenerife –38038–, otro supera los 15 mil habitantes, Santa Cruz de la Palma –38037–, también en la provincia de Santa Cruz de Tenerife aunque en diferente isla, y otros 4 superan los 5.000. Ningún municipio, sin embargo, se encuentra por debajo de los 300m de altitud promedio, siendo el más bajo de ellos precisamente Santa Cruz de Tenerife –38038–. Claramente las Islas Canarias destacan en lo abrupto del terreno (Goerlich 2022), 10 municipios de este listado pertenecen a Santa Cruz de Tenerife y otros 3 a Las Palmas. El resto se distribuyen entre León, Huesca, Asturias, Cantabria y Lleida.

Tabla 2. Municipios con mayor índice de rugosidad

	Municipios					Rugosidad		Pendiente
Provincia	Código	Nombre	Población	Altitud	Superficie	Indice	Base 100	%	°
Santa Cruz de Tenerife	38033	San Andrés y Sauces	4,213	935	43	1.392	133.75	96.8	44.1
Santa Cruz de Tenerife	38037	Santa Cruz de la Palma	15,446	975	43	1.391	133.67	96.7	44.0
Santa Cruz de Tenerife	38027	El Paso	7,745	1,088	136	1.334	128.16	88.2	41.4
Santa Cruz de Tenerife	38049	Valle Gran Rey	4,643	599	32	1.324	127.21	86.7	40.9
León	24116	Posada de Valdeón	424	1,524	164	1.320	126.88	86.2	40.8
Santa Cruz de Tenerife	38010	Buenavista del Norte	4,766	492	66	1.317	126.56	85.7	40.6
Santa Cruz de Tenerife	38007	Barlovento	1,921	801	44	1.314	126.32	85.3	40.5
Huesca	22057	Bielsa	481	1,920	203	1.294	124.37	82.2	39.4
Asturias	33008	Cabrales	1,942	997	238	1.286	123.54	80.8	38.9
Cantabria	39088	Tresviso	59	1,029	16	1.285	123.46	80.6	38.9
Lleida	25063	Canejan	89	1,600	49	1.284	123.35	80.5	38.8
Santa Cruz de Tenerife	38021	Hermigua	1,767	629	39	1.280	123.01	79.9	38.6
Palmas, Las	35001	Agaete	5,695	452	44	1.280	122.97	79.8	38.6
Asturias	33003	Amieva	616	815	115	1.277	122.74	79.4	38.5
Huesca	22078	Canfranc	599	1,912	72	1.273	122.32	78.7	38.2
Huesca	22182	Plan	274	1,749	93	1.272	122.23	78.6	38.2
Huesca	22230	Torla-Ordesa	299	1,893	185	1.268	121.85	78.0	37.9
Palmas, Las	35020	La Aldea de San Nicolás	7,516	413	124	1.267	121.75	77.8	37.9
Asturias	33050	Ponga	574	1,004	206	1.267	121.75	77.8	37.9
Palmas, Las	35025	Tejeda	1,865	1,014	104	1.261	121.17	76.8	37.5
Santa Cruz de Tenerife	38036	San Sebastián de la Gomera	9,261	445	113	1.261	121.16	76.8	37.5
Santa Cruz de Tenerife	38038	Santa Cruz de Tenerife	208,563	345	150	1.260	121.05	76.6	37.4
Santa Cruz de Tenerife	38016	Garafía	1,813	951	103	1.258	120.93	76.4	37.4
Huesca	22122	Hoz de Jaca	68	1,607	12	1.256	120.67	75.9	37.2
Santa Cruz de Tenerife	38050	Vallehermoso	2,917	593	109	1.255	120.65	75.9	37.2
Fuente: MDT05 del IGN, Padrón 2021 del INE y elaboración propia

A efectos de examinar la heterogeneidad del terreno es conveniente analizar intervalos al igual que hemos hecho con la altitud –Tabla 1–, aunque para la rugosidad resulta menos intuitivo fijar cortes adecuados. La literatura ha propuesto varios intervalos para el índice TRI de Riley, DeGloria y Elliot (1999), ya que este es el índice más popular. Sin embargo, TRI depende de las unidades de medida y de la resolución del MDT de partida, lo que dificulta la generalización de intervalos comparables (Goerlich y Cantarino 2010a). Nuestros índices, $TRI_1$ o $TRI_2$, no están exentos de esta limitación aunque carezcan de unidades de medida. Los resultados de Goerlich (2022) muestran claramente como estos índices son bastante sensibles a la resolución del MDT de origen, ya que la rugosidad del paisaje tiende a suavizarse rápidamente conforme disminuimos la resolución, y en consecuencia el valor de los índices disminuye rápidamente, piénsese que al reducir la resolución del MDT promediamos en dos dimensiones.³⁰

Puesto que existe una relación unívoca entre $TRI_1$ y la pendiente media del terreno, resulta más intuitivo efectuar cortes en términos de la pendiente en términos porcentuales.³¹ El análisis de los datos mostró que cuatro cortes de rugosidad con pendientes en el 10%, 20% y 35% eran adecuados para un análisis descriptivo. En términos de $TRI_1$ un valor de la pendiente del 10% representa un valor del índice de 1.005, el 20% un valor de 1.020 y el 35% un valor de 1.059. Aunque el valor central de corte, 20%, está bastante por debajo del valor nacional –el 28.8%–, la asimetría en la distribución –Figura 2– hace que el promedio simple de las pendientes a nivel municipal –24.1%³²– sea significativamente inferior al valor nacional, y que el valor mediano –21.4%–, que divide a los municipios en dos grupos iguales, esté relativamente cerca del punto de corte.³³

Con esta partición obtenemos los resultados de la Tabla 3. La columna $TRI_1$ es el valor del índice correspondiente a los municipios de cada grupo, y en consecuencia es consistente con el valor nacional que se ofrece en la última fila. Resulta interesante que, aunque la distribución porcentual de municipios por tramos no es muy diferente de la de la Tabla 1, la distribución de la población por tramos de rugosidad si es radicalmente diferente. Los grupos extremos concentran alrededor de un tercio de la población a partes iguales –16.2% en el primer tramo y 16.7% en el segundo–. Los otros dos cortes acogen, aproximadamente, a otro tercio de la población cada uno de ellos. No detectamos pues una preferencia por la planicie como si parece detectarse una preferencia por el valle en la Tabla 1, aunque ahora si observamos que la altitud tiende a crecer con la rugosidad. Claramente, la relación entre estas dos variables es globalmente no lineal y difícil de resumir mediante estadísticos simples.

Tabla 3. Indicadores municipales por cortes de rugosidad (Pendiente en %)

		Municipios1		Población		Superficie		Altitud
Rugosidad	$TRI_1$	n	%	N	%	km²	%	m	Base 100
Hasta 10%	1.003	1,491	18.2	7,662,719	16.2	72,008	14.2	621	90.81
(10%, 20%]	1.012	2,350	28.6	15,554,437	32.8	143,089	28.3	656	95.95
(20%, 35%]	1.036	2,523	30.7	16,266,095	34.3	177,049	35.0	686	100.28
Más de 35%	1.107	1,848	22.5	7,901,856	16.7	114,094	22.5	767	112.09
España	1.041	8,212	100.0	47,385,107	100.0	506,240	100.0	684	100.00
1 El número de municipios incluye los territorios mancomunados
Fuente: MDT05 del IGN, Padrón 2021 del INE y elaboración propia

Finalmente resulta interesante combinar ambas variables. El Mapa 3 resume la información de los Mapas 1 y 2 mediante un mapa bivariante³⁴ y no puede ser más claro. ¡Ciertamente los extremos se tocan! La cadenas montañosas identificadas en el Mapa 1 presentan elevada altitud y rugosidad –especialmente evidente en los Pirineos, Picos de Europa y la Isla de Tenerife– al tiempo que las depresiones presentan valores bajos de altitud y rugosidad –sobre todo la del Guadalquivir y en menor medida la del Ebro–, pero en una gran parte del territorio la relación entre ambas variables no es evidente. Así, podemos observar amplias zonas de la península con altitudes elevadas y rugosidad baja –la meseta castellana y amplias zonas de Castilla-La Mancha– al tiempo que muchas zonas costera de baja altitud presentan una rugosidad notable –la cornisa cantábrica y el País Vasco, pero también la mayor parte de la costa mediterránea con pocas excepciones, desde Andalucía hasta Cataluña–.

---

---

La Tabla 4 pone algunos números a las impresiones visuales del Mapa 3 y representa, con los cortes de las Tablas 1 y 3, la distribución bivariante de altitud y rugosidad en términos de municipios, panel (a), y población, panel (b).

El panel (a), que ofrece la distribución de municipios, indica claramente que los más altos no son necesariamente los más abruptos, ya que la mayoría de los municipios de más de 1000m de altitud tienen una rugosidad en el entorno del promedio nacional –el 38.4%–. Por otra parte, algo más del 11% de municipios presentan altitudes por encima del promedio, pero podrían ser considerados planos, ya que su pendiente media no supera el 10%, mientras que los municipios con menor altitud no necesariamente tienen una menor rugosidad.

El panel (b), que ofrece la distribución de la población, muestra, como hemos observado anteriormente, que la preferencia por altitudes bajas no se traduce directamente en una preferencia por el terreno llano. Casi la mitad de la población que vive en municipios con altitudes hasta los 300m lo hace en municipios con una rugosidad en el entorno del promedio nacional–esta celda es la moda de la distribución bivariante con casi la cuarta parte del total–, y un 15% lo hace en municipios con un elevado grado de rugosidad, y es que una parte importante de nuestras costas es relativamente abrupta. En el otro extremo de la distribución si observamos que la población que reside en municipios elevados tiende a hacerlo, con porcentajes crecientes, en municipios con elevada rugosidad, aunque las cifras absolutas son muy reducidas. La celda de mayor altitud/mayor rugosidad englobaría lo que podríamos realmente definir como poblaciones de montaña, unos 500 municipios que albergan alrededor del medio millón de personas, ¡pero el 40% de esta población se debe a Santa Cruz de Tenerife –38038–!

Tabla 4
(a) Distribución porcentual de municipios

	Rugosidad
Altitud	Hasta 10%	(10%, 20%]	(20%, 35%]	Más de 35%	España
Hasta 300m	3.91	5.15	6.30	3.45	18.80
(300, 700]	3.09	8.56	10.06	7.23	28.95
(700, 1000]	10.69	10.38	7.42	5.66	34.15
Más de 1000m	0.46	4.53	6.95	6.16	18.11
España	18.16	28.62	30.72	22.50	100.00

(b) Distribución porcentual de población

	Rugosidad
Altitud	Hasta 10%	(10%, 20%]	(20%, 35%]	Más de 35%	España
Hasta 300m	10.73	11.63	23.04	8.08	53.48
(300, 700]	2.99	16.85	8.06	5.99	33.89
(700, 1000]	2.44	3.82	2.28	1.49	10.02
Más de 1000m	0.01	0.53	0.96	1.11	2.61
España	16.17	32.83	34.33	16.68	100.00
Fuente: MDT05 del IGN, Padrón 2021 del INE y elaboración propia

6. Altitud, rugosidad y distribución de la población.

Es conocido que el proceso de concentración de la población que España viene experimentado desde hace más de un siglo ha producido un desplazamiento importante “de la montaña al valle” (Reher 1986, Goerlich y Mas 2009, Collantes 2004, Goerlich y Molla 2021). Cuantitativamente, “la montaña” se ha medido casi siempre a través de la altitud de los núcleos o del término municipal, y no se han tenido en cuenta otros factores derivados de una orografía adversa que también podrían identificar ciertas zonas como áreas de montaña.

Así pues, este último epígrafe examina la evolución histórica de la población, de 1900 a 2021, según la altitud, pero también de acuerdo con la rugosidad, a partir de las poblaciones censales homogéneas de Goerlich, Ruiz, Chorén y Albert (2015) a las que añadimos la población del Padrón de 2021.³⁵

La Figura 3 muestra la evolución histórica de la población por cortes altimétricos. Esta evolución es bien conocida (Goerlich y Mas 2009, Goerlich y Molla 2021) y muestra, fundamentalmente, la concentración de la población en la franja costera. En 1900 el 35.9% de la población se localizaba en municipios con una altitud media que no sobrepasaba los 300m, pero en 2021 era más de la mitad de la población –53.5%–. El crecimiento fue especialmente acelerado entre los censos de 1950 y 1981, 10 puntos porcentuales, coincidiendo en gran medida con la industrialización de España. El resto de franjas altimétricas pierden población. Esta pérdida es muy escasa en el intervalo entre los 300 y los 700m, a lo que sin duda contribuye el hecho de que Madrid pertenece al mismo –con una altitud media de 656m– y ha multiplicado su población por 6 –frente al 2.5 del promedio nacional–, alcanzando en 2021 los 3.3 millones de residentes. En términos relativos el intervalo entre los 700 y los 1000m de altitud pierde la mitad de su población, 10 puntos porcentuales, del 20.3% al 10.0%, y los municipios por encima de los 1000m de altitud promedio pasan de representar el 8.7% de la población a tan solo el 2.6%. Esta es la única franja que ve caer la población en términos absolutos, pasando de 1.6 millones de personas en 1900, hasta 1.2 millones en 2021. La mayor población de este intervalo fue alcanzada en 1950 con 1.9 millones de personas.

---

---

La evolución mostrada en la Figura 3 puede contraponerse a la que se observa en la Figura 4, y que es mucho menos conocida. Esta imagen recoge la distribución de la población según los cortes de rugosidad examinados en el apartado anterior en base a pendientes medias municipales derivadas del índice $TRI_1$ –Tabla 3–. La evolución es radicalmente diferente a la observada en la Figura 3. Aquí observamos una mayor estabilidad en la distribución de la población según los cortes que hemos hecho. En términos absolutos ninguno de los intervalos considerados ve disminuir sus efectivos demográficos, y en términos relativos las pérdidas de población se concentran mayoritariamente en los municipios con pendiente media superior al 35% –índice $TRI_1$ superior a 1.059–, que pasan de representar el 21.5% de la población en 1900 al 16.7% en 2021, aunque en términos absolutos experimentan un crecimiento de casi 3 millones de efectivos demográficos. El trasvase relativo de población se produce a zonas de rugosidad intermedia-baja, con pendientes entre el 10% y el 20%, que pasan de representar el 27.4% de la población en 1900 al 32.8% en 2021. Las otras dos franjas experimentan escasas variaciones, tan solo de 1.5 puntos porcentuales, de incremento en el caso de las áreas más planas y de disminución en el caso de los municipios con pendientes entre el 20% y el 35%. El resumen parece evidente, la rugosidad del terreno parece jugar un papel mucho más moderado en el proceso de concentración de la población que la altitud y una parte importante de ello se debe a que muchas zonas de baja altitud experimentan valores de rugosidad moderados o incluso elevados –Mapa 3–.

Así pues, ambos aspectos, altitud y rugosidad, parecen relevantes a la hora de catalogar lo que podríamos denominar como áreas de montaña, ya que ambas variables tienen contenido informativo diferente. De hecho, así lo entiende la Ley 25/1982, de agricultura de montaña, si bien, hasta donde nosotros conocemos, no existe una tipología oficial de estas áreas (ATECMA 2006), cuya cuantificación no resulta sencilla y cuyos umbrales parecen ser extremadamente sensibles a la resolución de la información de partida (Goerlich 2022).

---

---

Finalmente, a efectos comparativos, la Tabla 5 muestra la distribución bivariante de altitud y rugosidad en términos de municipios para la población en 1900, y que puede compararse con la existente en 2021 –Tabla 4(b)–. De ambas tablas se deduce que los movimientos de población a zonas de baja altitud, hasta los 300m, se produce independientemente de la rugosidad, puesto que en todos los cortes se observan incrementos relativos de población, siendo los más intensos en el intervalo alrededor del promedio nacional, 7 puntos porcentuales de incremento. El fenómeno contrario se observa en las zonas de altitud por encima de los 700m, donde observamos disminuciones de la población, en términos relativos, en todos los intervalos de rugosidad. Por el contrario, en el rango de altitudes intermedias, por encima de los 300 y hasta los 700m, que muestra una cierta estabilidad en la población –con una ligera caída de 1.2 puntos porcentuales–, se observan crecimientos en las zonas de baja rugosidad –hasta el 20%, pero relativamente importante en el rango entre el 10% y el 20%– y disminuciones en las zonas de rugosidad media y elevada –por encima del 20%, pero relativamente importante en el intervalo medio, entre el 20% y el 35%–.

Tabla 5. Distribución porcentual de población en 1900

	Rugosidad
Altitud	Hasta 10%	(10%, 20%]	(20%, 35%]	Más de 35%	España
Hasta 300m	6.70	7.71	16.00	5.48	35.89
(300, 700]	2.58	12.68	12.22	7.60	35.08
(700, 1000]	5.11	5.45	5.28	4.47	20.32
Más de 1000m	0.11	1.52	3.17	3.92	8.71
España	14.50	27.36	36.67	21.47	100.00
Fuente: MDT05 del IGN, Goerlich et al (2015) y elaboración propia

7. Reflexiones finales

Este trabajo ha ofrecido una actualización de los cálculos realizados en Goerlich y Cantarino (2010a). ¡Lo hemos actualizado todo, 👍! La información original, aumentando notablemente la resolución, los procesos de cálculo y el índice. El nuevo índice se basa en una estimación de la superficie del paisaje –superficie 3D– a nivel de celda de 5m $\times$ 5m –25m² de resolución–, es mucho más intuitivo, permite una transformación ordinalmente equivalente que acota el índice al intervalo [0, 1] si así se desea, puede transformarse a un indicador de la pendiente (media) del área en cuestión y es consistente en la agregación a partir de la agregación de sus componentes –las superficies 2D y 3D– y no a partir de la promediación de los índices elementales.

Las conclusiones sobre la relación entre altitud y rugosidad, y la relación entre estas variables y la distribución de la población no se alteran en lo sustancial respeto a lo que ya observamos en Goerlich y Cantarino (2010a) para el siglo XX, y que ahora hemos extendido hasta 2021. Altitud y rugosidad ofrecen información diferente sobre la orografía del paisaje, presentan una correlación relativamente baja, inexistente a nivel provincial, y los patrones de asentamiento de la población son diferentes respecto a ambas variables. Si los movimientos demográficos muestran una clara preferencia por las zonas de baja altitud esto es mucho menos evidente en relación a las zonas de baja rugosidad. Este resultado se debe, en gran parte, a que España tiene extensas zonas de baja altitud, próximas a la costa, con rugosidades elevadas –Mapa 3–.

Aunque esta aplicación hace uso del MDT05 de rugosidades derivado en Goerlich (2022), con tamaño de pixel de 25m², agregado a nivel de municipio. La elevada resolución de esta información permite su agregación a prácticamente cualquier escala geográfica y permite la definición de las denominadas áreas de montaña independiente de los lindes administrativos a partir de la combinación de información sobre altitudes y rugosidades.

---

Referencias

• Abramowitz, M. & Stegun, I. A. (1972) Handbook of mathematiccal functions with formulas, graphs and mathematical tables Dover Publications, New York, New York, USA.

• Aguilar Torres, F. J.; Aguilar Torres, M. A.; Agüera Vega, F.; Carvajal Ramírez, F. & Sánchez Salmerón, P. L. (2002) “Efectos de la morfología del terreno, densidad muestral y métodos interpolación en la calidad de los modelos digitales de elevaciones”. XIV Congreso Internacional de Ingeniería Gráfica. Santander. 5-7 de Junio de 2002.

• ATECMA (2006) Fundamentos para la definición de Zona de Montaña en España y estudio de la aplicación de la indemnización compensatoria en zonas de montaña Dirección General para la Biodiversidad. Ministerio de Medio Ambiente.

• Beasom, S. L.; Wiggers, E. P. & Giordono R. J. (1983) “A technique for assessing land surface ruggedness”. Journal of Wildlife Management. 47, 1163-1166.

• Collantes, F. (2004) “El declive demográfico de la montaña española (1850-2000) ¿Un drama rural?” Ministerio de Agricultura, Pesca y Alimentación. Madrid.

• Crawford, M. A. (2008) An Analysis of Terrain Roughness: Generating a GIS Application for Prescribed Burning. Master of Science Thesis. Texas Tech University.

• Evans, J. S. (2021) spatialEco. R package version 1.3-6.

• Evans, J. S.; Oakleaf J. & Cushman S. A. (2014) An ArcGIS Toolbox for Surface Gradient and Geomorphometric Modeling. Version 2.0-0 Consultado: 20-09-2022.

• Felicísimo, A. M. (1994) Modelos Digitales del Terreno. Introducción y aplicaciones en las ciencias ambientales Pentalfa, Oviedo 1994. Biblioteca de Historia Natural.

• GDAL/OGR contributors (2022) GDAL/OGR Geospatial Data Abstraction software Library.

• Gionfriddo, J. P. y Krausman, P. R. (1986) “Summer habitat use by mountain sheep”. Journal of Wildlife Management. 50, 331-336.

• Goerlich, F. J. (2022) “Superficie planimétrica versus superficie del paisaje en España. –superficie 2D versus superficie 3D–” Documento de Trabajo del Ivie. WP-2022-07.

• Goerlich, F. J. & Cantarino, I. (2010a) “Un índice de rugosidad del terreno a escala municipal a partir de modelos de elevación digital de acceso público”. Documentos de Trabajo 7-2010 (Julio). Fundación BBVA. Bilbao.

• Goerlich, F. J. & Cantarino, I. (2010b) “Rugosidad del terreno. Una característica del paisaje poco estudiada”. Documentos de Trabajo 10-2010 (Octubre). Fundación BBVA. Bilbao.

• Goerlich, F. J. & Pérez, P. (2021) LAU2boundaries4spain: R package providing LAU2 (municipalities) data geometries for Spain for 2002-2021 ROpenSpain. Repositorio en GitHub.

• Goerlich, F. J. & Mas, M. (2009) “Drivers of Agglomeration: Geography versus History”. The Open Urban Studies Journal, 2, 28-42.

• Goerlich, F. J. & Mollá, S. (2021) “Desequilibrios demográficos en España: evolución histórica y situación actual”. Presupuesto y Gasto Público 102, 31–54. Instituto de Estudios Fiscales.

• Goerlich, F. J.; Ruiz, F.; Chorén, P. & Albert, C. (2015) Cambios en la estructura y localización de la población: una visión de largo plazo (1842-2011) Primera edición. Fundación BBVA. Bilbao. España.

• Hijmans R. (2022a) terra: Spatial Data Analysis. R package version 1.6-17.

• Hijmans R. (2022b) raster: Geographic Data Analysis and Modeling. R package version 3.6-3.

• Hobson, R. D. (1972) “Surface roughness in topography: Quantitative approach”. Chapter 8 in R. J. Chorley (Ed.) Spatial Analysis in Geomorphology. Harper & Row, New York, USA. 221-245.

• Jennes, J. S. (2004) “Calculating landscape surface area from digital elevation models”. Wildlife Society Bulletin. 32, 3, (September), 829-839.

• Pebesma, E. (2018) “Simple Features for R: Standardized Support for Spatial Vector Data”. The R Journal. 10, 1, 439-446. doi: 10.32614/RJ-2018-009.

• Polidori, L. J.; Chorowicz, J. y Guillande, R. (1991) “Description of terrain as a fractal surface, and application to digital elevation model quality assessment”. Photogrammetric Engineering & Remote Sensing. 57, 1329-1332.

• R Core Team (2022) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria.

• Reher, D. S. (1986) “Desarrollo urbano y evolución de la población: España 1787-1930”. Revista de Historia Económica 4, 1, 39-66.

• Riley, S. J., DeGloria, S.D. & Elliot, R. (1999) “A Terrain Ruggedness Index that Quantifies Topographic Heterogeneity” & “Erratum”. Intermountain Journal of Science. 5, 1-4 (December), 23-27.

• Sappington, J. M.; Longshore, K. M. & Thompson, D. B. (2007) “Quantifying Landscape Ruggedness for Animal Habitat Analysis: A Case Study Using Bighorn Sheep in the Mojave Desert”. Wildlife Society Bulletin. 71, 1419-1426. doi: 10.2193/2005-723

• Wiggers, E. P. y Beasom, S. L. (1986) “Characterization of sympatric or adjacent habitats of 2 deer species in West Texas”. Journal of Wildlife Management. 50, 129-134.

• Wilson, M. F. J.; O’Connell, B.; Brown, C.; Guinan, J.C. & Grehan, A.J. (2007) “Multiscale Terrain Analysis of Multibeam Bathymetry Data for Habitat Mapping on the Continental Slope”. Marine Geodesy. 30, 3-35. doi: 10.1080/01490410701295962.

• Wickham et al., (2019) “Welcome to the tidyverse” Journal of Open Source Software. 4, 43, 1686. doi: 10.21105/joss.01686.