Google Layers en QGIS

Como casi siempre, trasteando con QGIS, me encuentro con una extensión bastante útil: Google Layers (elaborada por Ivan Mincik y Nungyao Lin). Para acceder a ella, Complementos -> Obtener Complementos de Python, añadir los repositorios de terceros y ubicar e instalar Google Layers (pueden filtrar la búsqueda con Google). La utilidad radica en que se puede utilizar la imaginería de Google Earth/Maps en los proyectos sin necesidad de georreferenciar previamente lo que reduciría los errores por propagación en los cambios de proyección. Los diferentes tipos de mapas a los que se puede tener acceso son Road Map, Mobile Map, Terrain Map, Satellite Map e Hybrid. Las condiciones de uso establecen que un mismo usuario no puede tener a acceso a más de 1000 imágenes por día. Los proyectos deberían configurarse con la proyección de Google (Google Mercator) cuyo código EPSG (European Petroleum Survey Group) es 900913. Los usuarios europeos tienen datos de muestra de Europa en la ruta del plugin correspondiente a .qgis/python/plugins/GoogleLayers/shp.

Para facilitar el uso del plugín y no malgastar el cupo de 1000 imágenes diarias en intentos de ensayo y error recomiendo crear un shape de puntos delimitador de nuestra zona de interés, con la proyección correspondiente, reproyectar a Google Mercator (Guardar como… en Propiedades de capa) y usarla como referencia para bajar las imágenes en un nuevo proyecto en Google Mercator. Una vez bajadas todas las imágenes necesarias, hacer un mosaico de ellas si es el caso, reproyectar nuevamente y cortar si se quiere que sea prolijo.

Para probar el uso del plugín usé un procedimiento diferente que me permitió, incluso, verificar que tan buena es la georreferenciación en QGIS. Reproyecté a Google Mercator una imagen de Google Earth que había georreferenciado en UTM WGS84 zona 19 N, la cargué en una vista de un Proyecto (unidades metros, EPSG: 900913) e invoqué el uso del plugín (Map Type: Satellite en Settings seguido de Complementos -> Google Layers -> Google Layers) en una vista con zum completo a la imagen. El resultado fue este:

En algunas zonas la concordancia es buena y en otros hay desplazamientos del orden de los 20 metros. Para bajar otra imagen nos desplazamos, por ejemplo, ligeramente hacia la izquierda y usamos nuevamente el plugín. Se obtuvo lo siguiente:

El empalme con la otra imagen de Google es exacto.

Zona buffer con QGIS

Supongamos que queremos visualizar los riesgos de afectación de obras civiles por las eventuales inundaciones de un río y que la zona de peligro está alrededor de 100 metros del eje del mismo. Esto puede ser establecido usando la opción Buffer en las herramientas de geoproceso. Para ejemplificar el procedimiento se va a cargar un ráster georreferenciado en UTM WGS 84, huso 19, correspondiente a un área cercana a la población de San Juan de los Morros, Estado Guárico. El río, fue digitalizado como polilínea y a su tabla atributiva se le añadió el campo buffer con el valor de 100 m; tal como se señala en la siguiente imagen:

Para determinar la zona buffer, Vectorial -> Herramientas de Geoproceso -> Buffer(s). En la ventana de diálogo que se presenta, marcar la casilla “Campo de distancia buffer” (que contiene el valor de 100), seleccionarlo en el menú desplegable más abajo y escoger el nombre del archivo de salida (también se puede escribir directamente la distancia buffer de 100 m digitándo el valor en la casilla correspondiente sin necesidad de usar el proceso anterior). Después de OK se verá algo similar a lo de la imagen siguiente:

El archivo vectorial de tipo polígono correspondiente al buffer se ha hecho transparente en un 50 % para poder observar que en la parte central-superior existen riesgos de afectación de casas construidas muy cerca del eje del río.

Dos tipos de máscara, con marcadores, en QGIS

Supongamos que queremos individualizar la zona de las lagunas en la imagen georreferenciada a continuación (gabanas3.tif):

existen varias formas de hacerlo. Dos de ellas pueden ser las siguientes:

Máscara 1 (ráster): Para que quede prolijo, se puede hacer un corte provisional de la zona de interés (lagunas.tif) con clipper (Ráster -> Extraction -> Clipper) y luego un poligono de extensión de la capa (cuadro_lagunas.shp con Vectorial -> Herramientas de Investigación) para una extracción por máscara (Ráster -> Projections -> Combar (Reproyectar)).

El resultado es igual al anterior marcando “valores sin datos” igual a 0. En este caso, si se carga cualquier imagen entonces como lo que está alrededor del cuadro de las lagunas es transparente se verá resaltado.

Si no se marca tal opción, aparecerá un recuadro negro que se elimina volviendo a reproyectar con “valores sin datos” igual a 0.

Máscara 2 (vectorial): Se hace el polígono (cuadro_gabanas.shp) de extensión de la capa para gabanas3.tif y se usa el polígono de extensión de la capa para lagunas, cuadro_lagunas.shp, para insertar un anillo con autoensamblado en modo de edición. Luego, en propiedades de la capa, se hace de color blanco tanto el cuerpo del polígono como el borde. Si aparece un reborde negro, se pueden mover ligeramente los nodos del polígono, por pares (click en cada nodo con la tecla Ctrl presionada hasta que se vean de color azul), en modo de edición. El resultado es similar en aspecto a la penúltima imagen con zum a la capa. Si se quiere un zum particular como éste:

entonces en Ver -> Nuevo marcador se le da un nombre (por ejemplo Lagunas) y cada vez que se quiera invocar esa zona particular Ver -> Mostrar marcadores y click en Lagunas y luego en zum a.

Corrección de geometría con QGIS

Hace algún tiempo bajé este shapefile de Venezuela de una página que hoy está caida:

Venezuela.zip

y sólo hasta ahora se me ocurrió comprobar la geometría. Para ello, con la capa cargada y seleccionada en una vista: Vectorial -> Herramientas de geometría -> Comprobar validez de geometría produce 76 errores. Como no son muchos se puede intentar corregirlos manualmente en modo de edición. Por tanto, se configuran las opciones de autoensamblado, por ejemplo, a 20 píxeles sin intersecciones (Configuración -> Opciones de autoeensamblado) se conmuta edición y se hace click en el primer error de la ventana de “Comprobar validez de geometría”. Automáticamente, QGIS hace zum en el sitio del problema:

Allí se ve que es una intersección, se hace click en la herramienta de nodos, barra de edición, y se desplaza ligeramente para eliminarla. Se graba y se invoca nuevamente la herramienta de comprobación de eliminación de geometría. En esta nueva ejecución el resultado es 72. El proceso se ejecuta sucesivamente hasta eliminarlas todas. Las más fáciles de eliminar son las intersecciones. Si embargo, un par de ellas son como las que sigue:

donde es necesario recurrir a la barra de edición avanzada para eliminar partes. El resultado fue el siguiente:

La falta de empalme no da error de geometría. Si se configuró en las opciones de autoensamblado “evitar intersecciones”, al arrastrar el nodo por detrás de la falta de empalme ésta se elimina automáticamente.

Digitalización en QGIS de polígonos con anillos e islas: paso de multipartes a partes sencillas

En este artículo vamos a pasar un shapefile de puntos a polígono con el complemento Point2One. Para acceder a el, Complementos -> Obtener complementos de Python y en la pestaña de Repositorios, clickear en “Añadir repositorios de terceros”, ir a la pestaña Complementos y ubicarlo por filtro para su instalación. Una vez instalado, se carga el shapefile de puntos y se ejecuta Complementos -> Points2One -> Points2One, cuidando de hacerlo sólo con las opciones por defecto. Esta es un imagen del resultado:

Ahora se van a añadir anillos en el polígono. En modo de edición, se hace click en el icono que está señalado en la imagen de abajo y se digitalizan, por ejemplo, dos de tales estructuras. Es equivalente a digitalizar, por ejemplo, cuerpos de agua dentro de un predio.

Ahora se van a añadir “islas”. Para ello se hace click en el icono de selección de partes individuales (el color vira a amarillo) y se digitaliza, por ejemplo, una estructura similar a la que está a la derecha luego de hacer click en el icono de añadir parte (ver imagen de abajo).

Después de salvar los cambios, la tabla atributiva sólo señala un elemento por lo que estamos en presencia de un polígono multipartes.

Para pasar de multipartes a partes sencillas Vectorial -> Herramientas de Geometría -> Multipartes a partes sencillas, escogiendo un nombre apropiado para el shapefile de salida antes de OK. En la imagen de abajo se observa que ahora el shapefile de salida presenta dos elementos en su tabla atributiva.

Para borrar los anillos del shapefile creado, click en el icono señalado en la imagen de abajo y luego en cualquier punto de las estructuras de abajo. En la imagen siguiente se observa como desaparecen tales estructuras.

Para eliminar las “islas” se tiene que volver al shapefile original, hacer click en el icono señalado en la imagen de abajo y luego en cualquier punto de la estructura a eliminar. Acepta para eliminación incluso la forma original situada a la izquierda (imagen superior). Sin embargo, se va a proceder a eliminar la forma de la derecha; tal como se observa en la imagen siguiente:

Para finalizar, vamos a hacer algunos cálculos de geometría.

Polígono con anillos:

área = 438460.385437 m², perímetro = 3929.578772 m

Polígono sin anillos:

área = 485841.678162 m², perímetro = 2727.177602 m

Esto representa una diferencia de áreas de 47381.29273 m² a favor del polígono sin anillos, tal como era de esperar, y de 1202.40117 m a favor del polígono con anillos, lo que indica que en este último caso se está contabilizando al valor total el perímetro de los anillos.

Hacer join a tablas externas en QGIS

Desde hace algún tiempo, quería descubrir este rasgo en QGIS. Llegó coincidencialmente a través de un comentario de Miguel Luis (Málaga, España), quien fue el que encontró como asignar el tipo de campo en la tabla a enlazar. Ya había descubierto, por mi cuenta, como cargar un archivo de texto delimitado por tabulador (copy/paste de Excel) que, en la vista de QGIS, se desplegaba como el equivalente a un Event Theme de ArcGis. Sin embargo, ese no es el formato de capa vectorial de tipo CSV (comma separated values) que reconoce QGIS. Para ello hay que exportarlo en “Guardar como…” e importarlo posteriormente con un archivo asociado, *.csvt, que es el que contiene la descripción del tipo de dato (Integer, Real, Double, String) de cada uno de los campos del archivo *.csv.

Para probar el procedimiento, se va a unir la siguiente tabla arbitraria (tabla_join.txt; sin espacios entre los campos):

ID,a,b,c
1,5.01,2.33,2.33
2,5.01,2.33,2.33
3,5.01,2.33,2.33
4,5.01,2.33,2.33
5,5.01,2.33,2.33

a un archivo arbitrario de puntos (punto.shp) con proyección UTM WGS 84 (zona 19). La tabla anterior se guarda como tabla_join.csv y se carga en la vista con el tabla_join.csvt asociado que contiene:

Integer,Real,Real,Real

sin espacios entre ellos. En la vista de QGIS se ve así:

donde se puede observar que ambas tablas, la del *.csv y la del Event Theme, son idénticas. Cargamos el shapefile, hacemos zoom a la capa y desplegamos su tabla atributiva:

Luego, hacemos click con el botón derecho sobre la capa punto y en el menú contextual seleccionamos Propiedades y en la pestaña Uniones añadimos con el botón de + la capa tabla_join.csv con los valores por defecto allí expresados:

Aquí se refleja la unión en la tabla de Propiedades:

y aquí el shapefile punto.shp con su tabla atributiva mostrando la unión.

Para que quede definitivamente enlazada, Guardar como…, por ejemplo, punto_union.shp.

En la tabla anterior se puede eliminar el campo NULL (que es el común) con Table Manager.

Georreferenciar una proyección cilíndrica equidistante del mundo, centrada en el Pacífico, con FWTools

Recientemente, alguien preguntó en un foro como se podía georreferenciar una imagen del mundo pero centrada en el Pacífico. Cuando comencé a estudiar tópicos sobre teledetección, recordé este tutorial:

http://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial

que incluye una data para trabajar con el mismo. Entre ellas está una imagen, world.png, que corresponde a una proyección cilíndrica equidistante del mundo de resolución 2048×1024.

Esta imagen se puede cortar en dos trozos de exactamente 1024×1024 con Gimp e invertirlas de posición. Calzarán perfectamente dando la impresión de que ahora están centradas en el Pacífico. El Gimp tiene una herramienta que permite posicionarla exactamente en el punto de corte; tal como se puede ver en la siguiente imagen:

Luego, se crea una vista nueva de 2048×1024 y se adiciona cada porción, en su nueva posición, como capas. El movimiento de las capas se puede hacer pixel a pixel y se garantiza empalme perfecto.

Con FWTools (el link de bajada está en el tutorial), la instrucción:

gdal_translate -a_srs WGS84 -a_ullr -180 90 180 -90 world.png geoworld.tif

asigna la proyección a la primera imagen y esta instrucción:

gdal_translate -a_srs WGS84 -a_ullr 0 90 360 -90 world_c.png geoworld_c.tif

a la segunda.

Ambas deben empalmar al cargarlas en QGIS; tal como se ve en la imagen siguiente: