Módulo m.cogo de GRASS para convertir líneas de rumbo y distancia a coordenadas cartesianas en QGIS

En este artículo:

Convertir Lineas de Rumbo y distancia a Coordenadas Cartesianas

señalé como era posible usar una hoja de cálculo para convertir líneas de rumbo y distancia a coordenadas cartesianas. Ahora se va a señalar como hacerlo utilizando el módulo m.cogo de GRASS dentro de QGIS. Se va a utilizar el mismo ejemplo del link anterior y se debe construir un archivo que debe tener este formato (la primera línea es un comentario y puede ser eliminada, si lo desean) y guardarlo con un nombre sugestivo (en nuestro caso rumbos-se.dat):

#<rumbo> <distancia>
N 13:51:8 E 320.71
N 55:18:8 E 214.09
N 79:50:8 E 168.97
N 13:45:8 E 170
N 28:20:0 W 25
S 74:27:38 E 99.63
S 12:58:0 W 67
S 86:5:52 E 105
S 1:30:8 W 379.5
S 53:20:52 E 280.75
S 12:23:32 E 315.97
S 10:47:59 E 481.09
S 46:2:35 E 153.68
S 28:26:11 E 84.37
S 32:38:2 E 501.2
S 30:17:16 E 562.44
S 65:35:11 W 626.98
N 34:54:42 W 393.55
N 38:19:52 W 162.77
N 37:25:52 W 510.32
N 37:32:52 W 489.34
N 37:33:52 W 525.6
N 13:41:8 E 479.93

Para acceder a la herramienta es necesario abrir cualquier mapaset y luego procurar en las herramientas de GRASS el módulo m.cogo (ver imagen siguiente). En la ventana de diálogo del módulo colocar la ruta al archivo de entrada (en nuestro caso /home/zeito/rumbos.dat en Name to input file), establecer el par de coordenadas iniciales (20000,20000 para nuestro ejemplo) y la ruta al archivo de salida.

Después de click en Ejecutar, en la pestaña Salida se indica la finalización satisfactoria y en el archivo de salida tenemos lo siguiente:

20076.783907 20311.382620
20252.801451 20433.252846
20419.119520 20463.071650
20459.532526 20628.198234
20447.667516 20650.203269
20543.655666 20623.512222
20528.621928 20558.220670
20633.378500 20551.075002
20623.429635 20171.705432
20848.667979 20004.109936
20916.476008 19695.501588
21006.620994 19222.932578
21117.249326 19116.260582
21157.424867 19042.070133
21427.706479 18619.992782
21711.369400 18134.324067
21140.450504 17875.180217
20915.216777 18197.905132
20814.265999 18325.588386
20504.089874 18730.825691
20205.874932 19118.796675
19885.441114 19535.423040
19998.989292 20001.727250

que es idéntico al resultado que se encuentra aquí:

Convertir Lineas de Rumbo y distancia a Coordenadas Cartesianas

Cabe señalar que la ejecución en modo GUI a través de QGIS funciona con las opciones por defecto. Si disponemos de un archivo con etiquetas (primera columna) como éste:

# <etiqueta> <rumbo> <distancia>
1 N 13:51:8 E 320.71
2 N 55:18:8 E 214.09
3 N 79:50:8 E 168.97
4 N 13:45:8 E 170
5 N 28:20:0 W 25
6 S 74:27:38 E 99.63
7 S 12:58:0 W 67
8 S 86:5:52 E 105
9 S 1:30:8 W 379.5
10 S 53:20:52 E 280.75
11 S 12:23:32 E	315.97
12 S 10:47:59 E	481.09
13 S 46:2:35 E 153.68
14 S 28:26:11 E	84.37
15 S 32:38:2 E 501.2
16 S 30:17:16 E	562.44
17 S 65:35:11 W	626.98
18 N 34:54:42 W	393.55
19 N 38:19:52 W	162.77
20 N 37:25:52 W	510.32
21 N 37:32:52 W	489.34
22 N 37:33:52 W	525.6
23 N 13:41:8 E 479.93

para procesarlo es necesaria la cónsola de GRASS. Es el primero en la lista de módulos (ver imagen siguiente):

La ejecución de este comando en la cónsola:

m.cogo -l input=/home/zeito/rumbos.dat output=/home/zeito/salida3.dat coord=20000,20000

produce entonces una salida con etiquetas al final de cada coordenada:

20076.783907 20311.382620 1
20252.801451 20433.252846 2
20419.119520 20463.071650 3
20459.532526 20628.198234 4
20447.667516 20650.203269 5
20543.655666 20623.512222 6
20528.621928 20558.220670 7
20633.378500 20551.075002 8
20623.429635 20171.705432 9
20848.667979 20004.109936 10
20916.476008 19695.501588 11
21006.620994 19222.932578 12
21117.249326 19116.260582 13
21157.424867 19042.070133 14
21427.706479 18619.992782 15
21711.369400 18134.324067 16
21140.450504 17875.180217 17
20915.216777 18197.905132 18
20814.265999 18325.588386 19
20504.089874 18730.825691 20
20205.874932 19118.796675 21
19885.441114 19535.423040 22
19998.989292 20001.727250 23

Tabla de colores para ráster en Grass desde QGIS

Recientemente, tuve problemas para visualizar un ráster DEM (y otros tipos) desde Windows en el ambiente GRASS de QGIS. En Seven fue imposible la equalización de histogramas de DEMs y en XP, para ráster NDVI, se “colgaba” el programa e impedía acceder nuevamente al mapaset por lo que había que crear uno nuevo. Se ve que a la integración todavía le falta pulir algunos detalles en este sistema operativo. Eso me llevó a movilizarme a Linux (Debian) donde verdaderamente no tiene desperdicio. Después de instalarme la documentación de grass mediante synaptic (grass-doc), pude acceder a información verdaderamente valiosa. No obstante, “deambulando” por los módulos de GRASS me encontré el correspondiente a r.colors.table donde se tiene ya predeterminadas una serie de opciones que no están disponibles ni en ERDAS, ArcGis o el propio QGIS. Estas son las siguientes:

aspect: aspect oriented grey colors
aspectcolr: aspect oriented rainbow colors
bcyr: blue through cyan through yellow to red
bgyr: blue through green through yellow to red
byg: blue through yellow to green
byr: blue through yellow to red
celsius: blue to red for degree Celsius temperature
corine: EU Corine land cover colors
curvature: for terrain curvatures (from v.surf.rst and r.slope.aspect)
differences: differences oriented colors
elevation: maps relative ranges of raster values to elevation color ramp
etopo2: colors for ETOPO2 worldwide bathymetry/topography
evi: enhanced vegetative index colors
gdd: accumulated growing degree days
grey: grey scale
grey.eq: histogram-equalized grey scale
grey.log: histogram logarithmic transformed grey scale
grey1.0: grey scale for raster values between 0.0-1.0
grey255: grey scale for raster values between 0-255
gyr: green through yellow to red
haxby: relative colors for bathymetry or topography
ndvi: Normalized Difference Vegetation Index colors
population: color table covering human population classification breaks
precipitation: precipitation color table (0..2000mm)
precipitation_monthly: precipitation color table (0..1000mm)
rainbow: rainbow color table
ramp: color ramp
random: random color table
rstcurv: terrain curvature (from r.resamp.rst)
rules: create new color table based on user-specified rules read from stdin
ryb: red through yellow to blue
ryg: red through yellow to green
sepia: yellowish-brown through to white
slope: r.slope.aspect-type slope colors for raster values 0-90
srtm: color palette for Shuttle Radar Topography Mission elevation
terrain: global elevation color table covering -11000 to +8850m
wave: color wave 

Para probar algunas de ellas usé el DEM ya considerado, el ráster de “sombras” producido a partir del DEM con r.shaded.relief (usando las opciones por defecto del módulo) y un NDVI producido con la calculadora ráster de QGIS a partir de las bandas 3 y 4 de Landsat 7 para una región de Calabozo, Guárico.

Con Windows Seven el ráster DEM no se puede equalizar (área totalmente negra) y desaparece de la vista irreversiblemente. No obstante, los resultados que se obtienen con r.watershed son correctos. Con Windows XP se equaliza a una opción por defecto pero cuando se trata de usar las opciones de r.colors.table el programa se “cuelga” e inutiliza irreversiblemente el acceso al mapaset y a las herramientas de GRASS; siendo necesario crear uno nuevo. No obstante, el acceso a los ráster allí almacenados es posible.

En Linux Debian el DEM se despliega a la opción por defecto del XP (aspect oriented grey colors); tal como se observa en la imagen 1. En la imagen 2 se observa un equalizado en tonos de grises y en la 3 con una opción equivalente a la de pseudocolor por defecto de QGIS; donde predominan los tonos azules para las elevaciones bajas. En la imagen 4 se tiene el ráster de sombras obtenido con r.shaded.relief y en la 5 y 6 las imágenes fueron colorizadas con una paleta que refleja una escala relativa con base a la altura y donde predominan los tonos en verde para las más bajas. La imagen 7 combina las imágenes 5 y la 6 para lograr efectos tridimensionales y las dos últimas están referidas a un ráster NDVI (índice de vegetación de diferencias normalizadas); donde la 9 ha sido coloreada con la paleta especial para este tipo de ráster. La zona con vegetación vigorosa aparece en verde muy oscuro y las zonas de quema (y también las nubes) con tonos rosa.

Generar escurrimiento superficial (“ríos”) con GRASS en QGIS

En el artículo anterior se consideró la “delimitación” arbitraria de una cuenca hidrológica con GRASS en QGIS. En éste se va a usar el módulo r.watershed para generar un ráster de escurrimiento superficial que luego será vectorizado en dos pasos mediante el uso sucesivo de los módulos r.thin y r.to.vect.line de GRASS desde QGIS. Las áreas de captación que se consideraron arbitrariamente fueron las de 10.000 y 100.000 Has; respectivamente 111111 y 1111111 celdas.

El detalle de la imagen del ráster con 111111 celdas para el tamaño mínimo de la cuenca, aumentada por el uso de una escala de aproximadamente 1:50.000, es el siguiente:

Después de usar consecutivamente los módulos r.thin y r.to.vect.line de GRASS, guardar como archivo vectorial de ESRI, pasarlo de partes sencillas a multipartes (por campo value generado por GRASS) y finalmente editarlo, se obtuvo lo siguiente después de añadirle las capas de ríos y cuerpos de agua de la región:

Se observa que un área de captación mínima de 10.000 Has está muy cerca de producir un escurrimiento superficial que simula la hidrografía de la región. Cuando se asume que el área de captación mínimo es de 100.000 Has, el escurrimiento superficial queda reducido de manera extrema; tal como se evidencia en la siguiente imagen (en rojo):

Por tanto, para afinar la delimitación de la cuenca mínima es necesario realizar varias corridas con valores entre 10.000 y 100.000 Has (transformados a número de celdas).

Delimitar una cuenca hidrológica con GRASS en QGIS

En el artículo anterior, escribimos sobre cómo crear una localización y conjunto de mapas para GRASS en QGIS. El ejemplo que se tomó fue con el objetivo de delimitar una cuenca hidrólógica en esa área. Allí se consideró el uso de QGIS en entorno Linux pero, para variar, se va a usar Windows XP (se usó inicialmente Seven pero, a pesar de que las herramientas de GRASS funcionan adecuadamente, los DEM no pueden visualizarse por ensanchamiento de contraste y “desaparecen” de la vista). En el caso de XP, la localización se creará en el directorio C:\Document and Settings\nombre_usuario\grassdata y llevará el mismo nombre: tiznados. El DEM importado fue el mismo del artículo anterior y se renombró para el complemento de GRASS como N09W068_dem_canoaUTM19. Tiene una superficie de algo más de un millón de hectáreas (1.225.469,61586 Ha); obtenida de la tabla atributiva del polígono de extensión de la capa (Vectorial -> Herramientas de investigación -> Polígono de extensión de la capa).

La delimitación se hará con base en el DEM y usando la herramienta de GRASS r.watershed (Complementos -> GRASS -> Abrir herramientas de GRASS). Cuando se abre la ventana de diálogo de la herramienta aparece en uno de los campos el “tamaño mínimo para cada cuenca (número de celdas)”. El ráster es de 30×30 por lo que una celda ocupa una superficie de 900 m². Para considerar tamaños mínimos de 10.000, 100.000 y 1.000.000 Has para cada cuenca entonces los números de celdas serán aproximadamente de 111111, 1111111 y 11111111; respectivamente.

Cuando se ejecuta la herramienta (cuidando, en cada caso, de clickar también en el botón que circunscribe el primer cuadrado rojo de la imagen anterior), el resultado es “vacío” para una cuenca de un millón de hectáreas (lo cual era previsible porque el DEM refleja una superficie similar) y produce lo siguiente:

para las superfices de 10.000 y 100.000 Has de tamaño mínimo para cada cuenca; respectivamente. En cada situación se han sobrepuesto los mismos shapefiles que representan los ríos y cuerpos de agua que se encuentran en la región. En el segundo caso se observa que la zona de color azúl oscuro central es la que contiene el río Tiznados; aunque este es un afluente del Río Portuguesa y se une con él aproximadamente a unos 87 Km en línea recta, hacia el SW, desde el borde inferior de la cuenca que lo incluye en mayor extensión.

Para vectorizar el ráster como un shapefile tipo polígono se utiliza la herramienta r.to.vect.area. El resultado es el siguiente:

después de colorizar las áreas por un categorizado del campo value (generado por GRASS). Se puede observar que es posible pasar de partes sencillas a multiparte para eliminar los registros repetidos y luego agregar columnas de geometría para estimar el área y el perímetro. Esto lo hice en el ámbito de QGIS; no de GRASS.

En el próximo artículo se ejecutarán y analizarán las otras opciones de r.watershed.

Crear una localización y conjunto de mapas para GRASS en QGIS

La integración con GRASS ya viene incluida en QGIS. Para acceder a ella es necesario activarla en Complementos -> Administrar complementos; marcando el correspondiente a GRASS. Como estamos trabajando en Linux vamos a aprovechar el directorio /usr/local/share/datagrass; creado en:

GRASS en Debian

para incorporar un ráster DEM (ASTER 30×30 = 900 m²) que ha sido bajado gratuitamente del servidor de la NASA (con registro previo obligatorio). No obstante, ese ráster ha sido descargado con una proyección original en coordenadas geográficas WGS 84 (EPSG: 4326) y se ha transformado a UTM 19N PSAD 56 (EPSG: 24819) con Ráster -> Proyecciones -> Combar (Reproyectar).

Se prevé trabajar en un sector denominado tiznados cuyo mapaset será dispuesto por el usuario jose. Para crear la localización en GRASS, a través de QGIS, se tiene que Complementos -> GRASS -> Nuevo directorio de mapas; donde aparecerá la ruta (/usr/local/share/datagrass) a la base de datos existente. En la siguiente ventana marcamos “Crear nueva localización” y escribimos tiznados. Definimos la proyección como PSAD 56 UTM19N (EPSG: 24819) y en la siguiente ventana Venezuela, Establecer:

Después de click en “Siguiente” se estableció a jose como el usuario del mapaset antes de finalizar la creación del directorio. En cada caso particular son libres de asignar el nombre que más le parezca.

Una vez creada la localización y el mapaset se va a importar el ASTER DEM en PSAD 56/UTM 19N con Complementos -> GRASS -> Abrir herramientas de GRASS donde en la pestaña Lista de módulos filtramos y seleccionamos r.in.gdal. La imagen ASTER en su proyección original tiene este nombre ASTGTM_N09W068_dem.tif. Se renombró a ASTGTM_N09W068_dem_canoa.tif cuando se reproyectó a PSAD 56/UTM19N y a demN09W068_UTM19N_canoa en la pestaña Opciones del módulo r.in.gdal. Cuando se ejecuta el módulo la imagen se carga automáticamente a la vista de QGIS al hacer click en Ver Salida.

Este es el resultado:

usando la paleta de pseudocolor para las propiedades de capa, pestaña Estilo.

Si se quiere usar este DEM en cualquier proyecto de QGIS entonces Complementos -> GRASS -> Añadir capa ráster de GRASS. Por otra parte, si se desea desde GRASS usar sus herramientas entonces es necesario abrir primero el directorio donde está el mapaset. Para ello Complementos -> GRASS -> Abrir directorio de mapas. Aquí tenemos un esbozo de la gran cantidad de ellas que se pueden usar desde QGIS:

Partes sencillas a multipartes en un shape de catastro hipotético con QGIS

El usuario Digd (España), en el foro de Gabriel Ortiz, refiere que tiene un shapefile de catastro (polígono) cuya tabla atributiva tiene múltiples campos y uno de ellos corresponde a la referencia catastral. El problema que tiene es que existen viviendas que abarcan varios polígonos y, en consecuencia, presentan la misma referencia catastral por lo cual requiere unir en un sólo proceso todos aquellos polígonos que presenten éste mismo valor.

Para ello, se vislumbra que en QGIS el paso de partes sencillas a multiparte (Vectorial -> Herramientas de geometría -> Partes sencillas a multiparte) puede ser una alternativa válida para resolver el problema. Para probar el procedimiento se considera un shapefile de catastro hipotético. Este corresponde a una cuadrícula rectangular donde se asume que los polígonos internos eran los predios y a la tabla atributiva se le adicionó un campo denominado ref_cat. A ese campo se le copió, de manera arbitraria, el número de filas y se escogieron varios grupos, aleatoriamente, para asignarle la misma ref_cat; de tal manera de simular viviendas que abarcan varios polígonos. Al ejecutar la herramienta de partes sencillas a multipartes, con ref_cat como campo ID único, el resultado fue el esperado; tal como se señala en la siguiente imagen:

En la gráfica, la selección amarilla corresponde a los polígonos con ID igual a 9 en el shapefile multiparte; cuya tabla atributiva es la segunda en la imagen. Su ref_cat es igual a 10. No obstante, en la tabla atributiva del shapefile de partes sencillas (arriba) se observa que esos tres mismos cuadrados tienen diferente ID (9, 10, 11) pero la misma ref_cat.

Para los otros grupos se obtuvo el mismo comportamiento.

CAD-Tools: puntos de intersección en QGIS

Una de las herramientas que extraño en QGIS es la posibilidad de añadir de manera masiva coordenadas X,Y en la tabla atributiva de los shapefile de puntos; algo que realizo con SAGA GIS. Como mejor es algo que nada, si se tiene shapefiles tipo polígono y se quiere obtener simultáneamente shapefiles tipo polilínea y de puntos ello es posible con la herramienta “puntos de intersección” de la barra de CAD-Tools. Con estos últimos, si el número de puntos no es muy elevado, es posible construir un archivo de coordenadas X,Y (formato csv) empleando la extensión Numerical Vertex Edit y el block de notas. En la imagen siguiente se ejemplifica el uso de la herramienta:

Se puede observar que el uso de la herramienta “puntos de intersección”, el primer par en la barra de CAD-Tools, produce automáticamente dos archivos: CadLayer Points y CadLayer Lines. Estos archivos no son permanentes por lo cual es inprescindible, para preservar los cambios, usar la opción “Guardar como…” del administrador de capas (los cambios no se preservan salvando el proyecto).

El primer punto de CadLayer Points fue producto de seleccionar los dos segmentos que lo delimitan con el icono de la izquierda en la herramienta y luego haciendo click en el icono de la derecha se materializa el punto (apareciendo CadLayer Points y CadLayer Lines en la vista de capas). En modo de edición de CadLayer Points, se usa la extensión Numerical Vertex Edit para copiar el contenido de su browser al blog de notas. Así se hace hasta completar todos los puntos.

Finalmente, se hace un join a CadLayer Points del archivo vectorial *.csv que contiene las coordenadas X,Y (el *.csvt contendría Integer,Real,Real) y se guarda con otro nombre para preservar los cambios.