Cómo crear una capa de puntos desde un fichero *.txt con PyQGIS

En esta entrada de blog vamos a ver un ejemplo práctico de cómo importar en QGIS un fichero *.txt (el cual contiene campos con coordenadas X e Y) utilizando Python y cómo añadir la capa a nuestra base de datos.

Como ejemplo voy a utilizar un fichero *.txt con los terremotos más importantes desde el año 2150 a.C, del Centro Nacional de Datos Geofísicos de la NOAA, el cual se puede descargar del siguiente enlace.

Si no funcionase el enlace, probad a guardarlo como *.txt desde este enlace, haciendo click derecho sobre la página-Guardar Como: http://www.qgistutorials.com/downloads/signif.txt.

Si abrimos el fichero, observamos que entre los diferentes campos tiene un campo “longitude” y otro “latitude”. También nos debemos fijar en la estructura, que vemos que en este caso está separado por tabuladores.

Para evitar fallos en la creación de la capa vamos a eliminar los tabuladores y los vamos a reemplazar por puntos y comas. Abrimos el editor de Python en QGIS (Ctrl+Alt+P y click en el icono de “Show editor” arriba a la izquierda) e importamos las librerías que nos van a hacer falta. Seguidamente, definimos la variable “basepath” con el directorio de trabajo, en mi caso “Carpeta de usuario\Documents\entradas\csv_pyqgis”. Definimos con el método .path.join() dónde tenemos almacenado el fichero *.txt (en la ruta de trabajo), el cual se denomina “signif.txt”. Abrimos el fichero con la función open(). Como segundo argumento le definimos el modo de apertura del fichero, que en este caso primero será de lectura (‘r’).

/**
import os
from pathlib import Path

basepath = Path().home().joinpath("Documents","entradas","csv_pyqgis").as_posix()
txt = os.path.join(basepath,'signif.txt')
with open(txt, "r") as f:
    datos = f.read()
*/

Una vez abierto el fichero de texto se procede a la lectura con el método .read(). Después, con el método .replace() especificamos que reemplace los tabuladores (\t) por puntos y comas (;). Abriremos nuevamente el fichero con .write() para escribir de nuevo todo el texto con los caracteres reemplazados.

/**
datos = datos.replace('\t', ';')
with open(txt, "w") as f:
    f.write(datos)
*/

Ahora el fichero tendrá este aspecto.

Ahora nos falta en nuestro código la parte de convertir este fichero de texto en una capa Shapefile. Creamos una variable denominada “uri”, la cual va a contener la ruta donde tenemos el fichero *.txt, a la cual le vamos a concatenar:

–delimiter: especificamos que el separador son puntos y comas.

–crs: especificamos el EPSG del sistema de coordenadas.

–yField: especificamos el nombre del campo que contiene las coordenadas Y en el *.txt.

–xField: especificamos el nombre del campo que contiene las coordenadas X en el *.txt.

A continuación, creamos la capa con la clase QgsVectorLayer, donde especificamos en los parámetros la ruta (uri), el nombre que va a tener la capa (‘terremotos’) y que va a ser una capa delimitada por texto (‘delimitedText’). Añadimos la capa al proyecto con el método .addMapLayer.

/**
uri="file:///"+txt+"?delimiter=;&crs=EPSG:4326&yField=latitude&xField=longitude"
capa = QgsVectorLayer(uri, 'terremotos', 'delimitedtext')
QgsProject.instance().addMapLayer(capa)
*/

Si ejecutamos el código veremos que nos añade a nuestro proyecto una capa que se llama terremotos. Sin embargo, esta capa es temporal. Si cerramos el proyecto perdemos la capa, por lo que vamos a exportar la capa a formato Shapefile. Para ello, con el método .writeAsVectorFormat de la clase QgsVectorFileWriter exportaremos nuestra capa con extensión ESRI Shapefile a nuestra carpeta de trabajo, utilizando el mismo sistema de coordenadas de la capa temporal.

/**
QgsVectorFileWriter.writeAsVectorFormat(capa, basepath+"/terremotos.shp", "CP1250", capa.crs(), "ESRI Shapefile")
*/

Si cargamos nuestro fichero Shapefile veremos las entidades de punto que representan los terremotos.

El código completo del ejercicio en Python es el siguiente.

/**
import os
from pathlib import Path

basepath = Path().home().joinpath("Documents","entradas","csv_pyqgis").as_posix()
txt = os.path.join(basepath,'signif.txt')

with open(txt, "r") as f:
    datos = f.read()

datos = datos.replace('\t', ';')
with open(txt, "w") as f:
    f.write(datos)

uri="file:///"+txt+"?delimiter=;&crs=EPSG:4326&yField=latitude&xField=longitude"
capa = QgsVectorLayer(uri, 'terremotos', 'delimitedtext')
QgsProject.instance().addMapLayer(capa)
QgsVectorFileWriter.writeAsVectorFormat(capa, basepath+"/terremotos.shp", "", capa.crs(), "ESRI Shapefile")
*/

El siguiente objetivo de esta entrada es importar esta capa a nuestra base de datos de PostgreSQL. Para ello, lo primero que debemos hacer es conectarnos a la base de datos. Utilizaremos el “Postgis shapefile import/export manager”, que se nos instala junto con PostGIS. En Windows lo podemos encontrar pulsando la tecla Windows. El nombre puede variar en función de la versión de PostGIS que tengamos instalada. En el caso de PostGIS 3:

Una vez abierto hay que establecer la conexión a la base de datos. Para ello hacemos click en “View connection details…”. Tendremos que rellenar una serie de datos.

–Username: nombre de usuario de PostGIS.
–Password: nuestra contraseña de PostGIS.
–Server Host: el nombre del servidor donde está alojado nuestra base de datos. En mi caso “localhost” al ser un servidor local. A la derecha especificamos el puerto de conexión a la base de datos. Por defecto, el puerto es 5432, aunque en mi caso tengo establecido 5433.
–Database: nombre de la base de datos a la cual le vamos a añadir la capa. Muy importante que nuestra base de datos tenga la extensión de PostGIS. Si no la tenemos activada, debemos ir a nuestra base de datos y escribir lo siguiente en el editor de consultas:

CREATE EXTENSION postgis;

Si la conexión se ha podido establecer exitosamente, aparecerá el siguiente mensaje abajo.

Le damos a “Add File” y seleccionamos el Shapefile “terremotos.shp”. Es importante especificar el nombre de la columna que contiene la geometría (“geom”). Especificamos también el SRID. Una vez rellenados estos datos, haceos click en “Import”.

Nos aparecerá en la parte inferior un mensaje de que se ha completado la exportación.

Si vamos a PgAdmin veremos que tenemos una nueva tabla llamada terremotos con el campo “geom”. Podemos hacer una consulta básica dentro de nuestra tabla, como por ejemplo mostrar aquellos terremotos que hubo en el año (year) 1920 con una intensidad (intensity) superior a 7 (medida en la Escala sismológica de Mercalli, que va del 1 al 12). Se seleccionan los campos de intensity, location_n (localización) y geom.

Para terminar, si hacemos click en el icono de la vista del campo geom nos aparecen señalados en el mapa los terremotos de intensidad mayor a 7 que hubo en el año 1920.