¿Cómo reproyectar 500 archivos CSV de manera eficiente y fácil usando QGIS?

11

Lo sé, mi pregunta es similar a algunas anteriores en este sitio.

Tengo muchos archivos CSV (coordenadas geográficas) para importar a qgis (y luego convertirlos), y la forma habitual no es la mejor manera de hacerlo (demasiado tiempo).

Tengo casi 500 archivos CSV (coordenadas wgs84) y esto es lo que quiero hacer:

  1. Importe todos los archivos CSV a la vez a QGIS
  2. Proyectarlos
  3. Exportarlos a archivos CSV (nuevamente) pero con diferentes coordenadas (conversión a UTM33N)

Estoy tratando de entender cómo usar la consola de Python, pero no sigo adelante :(

¿Alguien puede explicarme cómo lograrlo paso a paso?

Raquel Ribeiro
fuente
mira mi respuesta a continuación. el problema ya estaba resuelto y explicado
Generic Wevers
2
¿Y por qué es ese duplicado con el marcado? Tal vez el OP intenta aprender pyqgis y cómo usar python si considera sus negrita.
nickves
Por favor especifique su pregunta. ¿Desea no cargarlos manualmente en QGIS? ¿Quieres convertirlos a otro formato? ¿Cuál es exactamente su pregunta?
bugmenot123
1. Importar todos los archivos en un proceso a qgis 2. proyectarlos 3. exportarlos nuevamente como csv pero en coordenadas utm
Raquel Ribeiro
cat * .csv> one_file.csv (o cualquiera que sea el equivalente de Windows) combinará todos sus archivos csv en uno. 500 realmente no es un número tan grande :-)
John Powell

Respuestas:

15

Si está buscando reproyectar archivos csv desde la Consola Python en QGIS, puede usar el siguiente script. Todo lo que necesitaría cambiar son las tres rutas que se mencionan en los comentarios.

Esencialmente, el script importa sus archivos csv en QGIS como archivos de forma (suponiendo que sus campos geométricos se denominen Xy Y). Luego usa los algoritmos qgis:reprojectlayery qgis:fieldcalculatorde Processing Toolbox para reproyectar y actualizar los campos Xy Ycon las nuevas coordenadas. Luego los guarda en una carpeta y los convierte en archivos csv en la ruta que especifique. Entonces, al final, ha actualizado los archivos de forma y los archivos csv en carpetas separadas.

import glob, os, processing

path_to_csv = "C:/Users/You/Desktop/Testing//"  # Change path to the directory of your csv files
shape_result = "C:/Users/You/Desktop/Testing/Shapefile results//"  # Change path to where you want the shapefiles saved

os.chdir(path_to_csv)  # Sets current directory to path of csv files
for fname in glob.glob("*.csv"):  # Finds each .csv file and applies following actions
        uri = "file:///" + path_to_csv + fname + "?delimiter=%s&crs=epsg:4326&xField=%s&yField=%s" % (",", "x", "y")
        name = fname.replace('.csv', '')
        lyr = QgsVectorLayer(uri, name, 'delimitedtext')
        QgsMapLayerRegistry.instance().addMapLayer(lyr)  # Imports csv files to QGIS canvas (assuming 'X' and 'Y' fields exist)

crs = 'EPSG:32633'  # Set crs
shapefiles = QgsMapLayerRegistry.instance().mapLayers().values()  # Identifies loaded layers before transforming and updating 'X' and 'Y' fields
for shapes in shapefiles:
        outputs_0 = processing.runalg("qgis:reprojectlayer", shapes, crs, None)
        outputs_1 = processing.runalg("qgis:fieldcalculator", outputs_0['OUTPUT'], 'X', 0, 10, 10, False, '$x', None)
        outputs_2 = processing.runalg("qgis:fieldcalculator", outputs_1['OUTPUT_LAYER'], 'Y', 0, 10, 10, False, '$y', shape_result + shapes.name())

os.chdir(shape_result)  # Sets current directory to path of new shapefiles
for layer in glob.glob("*.shp"):  # Finds each .shp file and applies following actions
        new_layer = QgsVectorLayer(layer, os.path.basename(layer), "ogr")
        new_name = layer.replace('.shp', '')
        csvpath = "C:/Users/You/Desktop/Testing/CSV results/" + new_name + ".csv"  # Change path to where you want the csv(s) saved
        QgsVectorFileWriter.writeAsVectorFormat(new_layer, csvpath, 'utf-8', None, "CSV")   

¡Espero que esto ayude!

Joseph
fuente
2
Gran respuesta: ¡lo tienes todo allí! Una pregunta si no le importa: ¿todavía tiene que agregar / eliminar capas en QgsMapLayerRegistry incluso si hace cosas desde la consola de Python?
nickves
1
@nickves - Jaja muchas gracias amigo! Es posible que no tenga que agregar / quitar capas (estoy seguro de que el script se puede reducir drásticamente). No soy un experto, pero lo probaré más tarde y te responderé. A menos que pueda proporcionar un guión mucho más ordenado, en cuyo caso debe publicarlo como respuesta, lo votaría :)
Joseph
@nickves - ¡Gracias de nuevo por tu sugerencia amigo! El código ha sido editado para evitar agregar / eliminar las capas por segunda vez :)
Joseph
@RaquelRibeiro - ¡De nada! Me alegro de que haya sido útil :)
Joseph
@Joseph, ¿puedo preguntarte algo otra vez? En las líneas 14 y 15, los números: 0, 10, 10 definen exactamente qué. (las coordenadas de salida tienen demasiados ceros a la derecha y quiero minimizarlas)
Raquel Ribeiro
8

Una solución rápida para transformar un archivo separado por espacios que contiene "lon lat" en WGS84 a UTM33N pero no obtiene ningún otro dato:

#!/bin/bash
#
for i in $( ls *.csv ); do
    gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32633 < ${i} > utm${i}
done

Eso funciona y conserva el orden de los datos, así que ¿quizás otro bucle que utilice, por ejemplo, awk para combinar los datos descriptivos con las coordenadas?

Editar. Debido a los comentarios desordenados que hice a continuación, editaré la respuesta aquí.

El siguiente script debería hacer el trabajo de leer múltiples archivos csv, agregando nuevas columnas de coordenadas a cada archivo.

#!/bin/bash
#
for i in $( ls *.csv ); do
 paste -d',' ${i} <(awk -v OFS="," -F " " 'NR>1 {print $1 " " $2}' ${i} | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32633 | awk '{gsub(" ",",",$0); print $0}' | /usr/local/bin/sed "1i\X,Y,Z") > utm${i}
#
 #paste -d',' ${i} <(awk -v OFS="," -F " " 'NR>1 {print $1 " " $2}' ${i} | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32633 | awk '{gsub(" ",",",$0); print $0}' |sed "1i\X,Y,Z") > utm${i}
#
done

En OSX necesitará instalar la última versión (2009) de sed y usar la primera línea no comentada en el bucle. Para Linux, comente el primero y use el segundo. Ajústelo de -F " "acuerdo con el formato del separador en sus archivos csv, por ejemplo, -F ","para comas separadas. También tenga en cuenta que la transformación de elevación es al elipsoide, no al geoide, así que asegúrese de transformar las alturas en consecuencia.

mercergeoinfo
fuente
Acabo de recordar hacer algo similar hace un tiempo y publicar una solución en mi blog. Está escrito para Mac pero está basado en bash. La mayor diferencia es el problema con sed en OS X, que trato al final del post: mercergeoinfo.blogspot.se/2014/01/…
mercergeoinfo
El último comentario fue un poco desordenado. Use esta línea en el script bash anterior para recorrer todos los archivos. paste -d',' ${i} <(awk -v OFS="," -F " " 'NR>1 {print $1 " " $2}' ${i} | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32633 | awk '{gsub(" ",",",$0); print $0}' | /usr/local/bin/sed "1i\X,Y,Z") > utm${i}Reemplace / usr / local / sed con solo sed si no está en OSX. Esto no es ideal si sus archivos csv están separados por espacios, como se supone en la línea anterior, pero funciona. Si tiene una coma separada, cambie -F " "a-F ","
mercergeoinfo
Me pregunto por qué está el código actualizado en los comentarios y no actualizó su respuesta anterior. El código en el comentario es realmente difícil de leer. ¿Ves el enlace de edición debajo de tu respuesta?
Miro
Sí, pero en realidad no fue una actualización, más como una adicional. Bastante desordenado, estoy de acuerdo. Supongo que debería actualizar la respuesta original. Gracias
mercergeoinfo
7

Usar qgis o incluso OGR es excesivo para esto.
Uso pyproj( https://pypi.python.org/pypi/pyproj ) combinado con el escritor csv de python y algunos trucos de la biblioteca estándar. ¡No necesita instalar nada más que pyprojpara esto!

import csv
import pyproj
from functools import partial
from os import listdir, path

#Define some constants at the top
#Obviously this could be rewritten as a class with these as parameters

lon = 'lon' #name of longitude field in original files
lat = 'lat' #name of latitude field in original files
f_x = 'x' #name of new x value field in new projected files
f_y = 'y' #name of new y value field in new projected files
in_path = u'D:\\Scripts\\csvtest\\input' #input directory
out_path = u'D:\\Scripts\\csvtest\\output' #output directory
input_projection = 'epsg:4326' #WGS84
output_projecton = 'epsg:32633' #UTM33N

#Get CSVs to reproject from input path
files= [f for f in listdir(in_path) if f.endswith('.csv')]

#Define partial function for use later when reprojecting
project = partial(
    pyproj.transform,
    pyproj.Proj(init=input_projection),
    pyproj.Proj(init=output_projecton))

for csvfile in files:
    #open a writer, appending '_project' onto the base name
    with open(path.join(out_path, csvfile.replace('.csv','_project.csv')), 'wb') as w:
        #open the reader
        with open(path.join( in_path, csvfile), 'rb') as r:
            reader = csv.DictReader(r)
            #Create new fieldnames list from reader
            # replacing lon and lat fields with x and y fields
            fn = [x for x in reader.fieldnames]
            fn[fn.index(lon)] = f_x
            fn[fn.index(lat)] = f_y
            writer = csv.DictWriter(w, fieldnames=fn)
            #Write the output
            writer.writeheader()
            for row in reader:
                x,y = (float(row[lon]), float(row[lat]))
                try:
                    #Add x,y keys and remove lon, lat keys
                    row[f_x], row[f_y] = project(x, y)
                    row.pop(lon, None)
                    row.pop(lat, None)
                    writer.writerow(row)
                except Exception as e:
                    #If coordinates are out of bounds, skip row and print the error
                    print e
castillo-blord
fuente
Me doy cuenta de que el póster no tiene mucha experiencia con Python. No uso QGIS regularmente, entonces, ¿alguien con más experiencia en esa plataforma podría explicar dónde está instalado Python? El póster debe hacer de este un script independiente y probablemente ejecutarlo desde IDLE. No tengo una instalación actual, por lo que no sé si pyprojdebe instalarse por separado para el póster o si ya está allí.
blord-castillo
1
Nunca utilicé la función parcial antes. Lo haré de ahora en adelante. +1
nickves
4

No necesitas python. Simplemente use la línea de comando y ogr2ogr. En su caso, lo más importante es el parámetro -t_srs srs_def.

Esto ya se explica en esta respuesta a ¿Cómo puedo convertir un archivo de Excel con columnas x, y en un archivo de forma?

ACTUALIZACIÓN No tengo tiempo para escribirle su código completo. Pero el problema será que necesita un poco más de código en Python de lo que piensas.

Su principal problema será que trabajar con archivos csv no es tan cómodo como usar shapefiles. Por lo tanto, primero deberá convertir el csv a la forma que necesita el archivo VRT. Esto se explica en el primer enlace. Aquí deberá escribir un script de Python que recorra sus archivos y que genere automáticamente los archivos vrt.

Este es un guión que utilicé yo mismo. Tienes que probar si te funciona. Ya incluí la conversión de WGS 84 a UTM 33N

from os import listdir, stat, mkdir, system
path = "your path here"
out_path = "your output path here"
files = filter(listdir(path), '*.csv') #for Python 3.x
# files= [f for f in listdir(path) if f.endswith('.csv')] #for Python 2.7

for x in range(len(files)):
    name = files[x].replace('.csv', '')
    # 2. create vrt file for reading csv
    outfile_path1 = out_path + name + '.vrt'
    text_file = open(outfile_path1, "w")
    text_file.write('<OGRVRTDataSource> \n')
    text_file.write('    <OGRVRTLayer name="' + str(name) + '"> \n')
    text_file.write('        <SrcDataSource relativeToVRT="1">' + name + '.csv</SrcDataSource> \n')
    text_file.write('        <GeometryType>wkbPoint</GeometryType> \n')
    text_file.write('        <LayerSRS>WGS84</LayerSRS> \n')
    text_file.write('        <GeometryField encoding="PointFromColumns" x="Lon" y="Lat"/> \n')
    text_file.write('        <Field name="Name" src="Name" type="String" /> \n')
    text_file.write('    </OGRVRTLayer> \n')
    text_file.write('</OGRVRTDataSource> \n')
    # 3. convert csv/vrt to point shapefile
    outfile_path2 = out_path + name + '.shp'
    command = ('ogr2ogr -f "ESRI Shapefile" -t_srs EPSG:32633' + outfile_path2 + ' ' +  outfile_path1)
    system(command)

Debe ajustar los parámetros para Nombre de campo , src , x e y de acuerdo con su archivo csv.

ACTUALIZACIÓN2

Después de pensar un poco, me pregunto ¿por qué quieres usar QGIS? Puede usar un script de Python como este para convertir directamente sus coordenadas de WGS a UTM. En este caso es un simple csv abierto, leer coordenadas, transformar coordenadas y guardarlo en un nuevo archivo.

Wevers genéricos
fuente
creo que esto no es lo que estoy buscando ... Tengo casi 500 archivos csv (coordenadas wgs84) y esto es lo que quiero hacer: 1. Importar todos los archivos csv de una vez a q gis 2. proyectarlos 3. exportarlos a archivos csv (nuevamente) pero con diferentes coordenadas (conversión a utm33N)
Raquel Ribeiro
Creo que necesito un proceso por lotes o algo así para hacerlo ...
Raquel Ribeiro
44
pero ¿por qué quieres hacer eso? 1. puede hacer lo mismo (lo que describió) desde la línea de comandos sin qgis. 2. puede hacer esto en modo por lotes. 3. en python es casi lo mismo. también usarías ogr2ogr
Generic Wevers
2
"Simplemente" usar la línea de comando realmente no es una respuesta. La línea de comando nunca es fácil de usar si no tiene idea de cómo hacerlo. Y realmente no puedo encontrar la solución en la respuesta vinculada. ¿Por qué no simplemente darle al pobre tipo un lote de ejemplo con ogr2ogr, y todo estaría bien?
Bernd V.
1
ok, 1. puedes leer gis.stackexchange.com/help/how-to-ask . después de eso y 5 minutos de google admitirás que la pregunta está muy poco investigada y puede resolverse con respuestas ya dadas. 2. Si aún no se puede resolver, supongo que todos estarán encantados de ayudar. pero como soy una buena persona, daré más pistas.
Wevers genéricos