Procesado de mapas vectoriales con R

Siempre me ha interesado todo lo que tenga que ver con mapas, proyecciones cartográficas y representaciones de información geográfica, en definitiva GIS, pero hasta ahora no había tocado el tema con R. Sacando partido de unos datos de la Comunidad de Madrid vamos a practicar un poco con mapas vectoriales.

Usaremos datos vectoriales en formato shapefile obtenidos de los Datos abiertos de la Comunidad de Madrid. En concreto se trata de la altimetría de la provincia en curvas de nivel separadas 20m, y consiste en una lista de 7.239 objetos de tipo SpatialLinesDataFrame que describen con gran precisión la topografía madrileña. Extrayendo los puntos de los que se componen las líneas mediante una conversión a SpatialPointsDataFrame, obtenemos 888.804 puntos lo que da una definición promedio de las curvas de nivel de 123 puntos/línea.

Las curvas de nivel (SpatialLinesDataFrame), pueden superponerse fácilmente a su descomposición en puntos individuales (SpatialPointsDataFrame). Aquí vemos la combinación de ambos datos en el área de 2km x 2km entorno a Manzanares el Real, lugar señalizado por las coordenadas UTM obtenidas de Google Maps (hacer clic para ver en formato PDF):

~~~

Representamos directamente el mapa completo con la función spplot() del paquete raster. Señalamos la localización de Manzanares sobre la que luego haremos zoom (hacer clic para ver imagen en formato PDF de alta definición):

Repetimos el dibujo pero esta vez con el mapa de colores relativo a la altitud del Puig Campana (1.406m), cima alicantina que usaremos como referencia con los datos raster en el siguiente artículo. Si todo va bien las regiones coloreadas deberían cuadrar perfectamente con las que obtendremos, y se comprueba que es así (clic para PDF):

Ahora realizamos un recorte cuadrado entorno a Manzanares el Real y representamos de nuevo. Podemos ver una amplia zona sin curvas de nivel correspondiente al embalse de Santillana, y a espaldas de Manzanares la zona de la Pedriza, con una gran densidad de curvas por lo escarpado del terreno (hacer clic para PDF en alta resolución):

Un archivo vectorial se puede rasterizar para obtener una interpolación de cotas asociada a la matriz del terreno. Con la función interpolate() del paquete raster obtenemos el siguiente modelo de elevación del terreno:

Al disponer del mapa de cotas en formato matricial ya podemos dibujarlo en 3D. La definición sin embargo no es tan buena como cuando partimos de un raster ya que la interpolación de las curvas de nivel resulta bastante pobre (notar las zonas de valor constante en forma de "bancales"):

~~~

A continuación vamos a calcular la distribución de altitudes en la Comunidad de Madrid partiendo de los datos vectoriales. La primera tentación sería acudir directamente a las cotas del fichero vectorial (recordemos que se trata de 7.239 curvas caracterizadas por su altitud), pero esto sería altamente erróneo al no estar garantizado que las altitudes del terreno estén uniformemente representadas en las curvas de nivel. De hecho no es así al existir mayor profusión de curvas en las zonas montañosas que en cotas más bajas, lo que arrojaría una distribución de altitudes desplazada a la derecha.

La solución correcta si queremos partir del fichero vectorial pasa por rasterizar éste, y una vez trasladadas las cotas vectoriales a una matriz espacial calcular la distribución sobre la misma.

Como la interpolación de valores vectoriales se lleva a cabo para toda la cuadrícula que contiene la provincia de Madrid, eliminar el efecto de la interpolación fuera de Madrid requiere un paso extra: hemos de enmascarar el fichero vectorial rasterizado con otro equivalente que proporcione la silueta provincial. Este fichero shapefile poligonal lo obtenemos del Centro Nacional de Información Geográfica:

Hecho el cruce de la interpolación contra el anterior perfil vectorial de la provincia con la función mask() del paquete raster, tenemos un raster con una buena definición de cotas y definido solo en el área de interés. Puede verse la interpolación en forma de regiones tipo Voronoi generada fuera de Madrid, al no disponerse de datos vectoriales en esas zonas:

Con el dato interpolado solo para la provincia de Madrid, resulta trivial calcular la distribución de altitudes en ella:

El cálculo de altitud media de 819m se separa solo 2m del que aparece en este ranking provincial, así que parece que nuestra rasterización funcionó bastante bien.

También podemos tener un 3D solo de Madrid a partir de los datos vectoriales que hemos rasterizado:

~~~

En el artículo 'Procesado de mapas raster con R' se tratan los mapas de tipo raster sobre la misma región geográfica.

Repositorio con el código R y archivos auxiliares (incluyendo los ficheros shapefile): GitHub.