Comment reprojeter 500 fichiers CSV efficacement et facilement à l'aide de QGIS?

11

Je sais, ma question est similaire à certaines anciennes sur ce site.

J'ai beaucoup de fichiers CSV (coordonnées géographiques) à importer dans qgis (puis à les convertir), et la manière habituelle n'est pas la meilleure façon de le faire (trop longue).

J'ai près de 500 fichiers CSV (coordonnées wgs84) et voici ce que je veux faire:

  1. Importez tous les fichiers CSV à la fois dans QGIS
  2. Projetez-les
  3. Exportez-les dans des fichiers CSV (encore) mais avec des coordonnées différentes (conversion en UTM33N)

J'essaie de comprendre comment utiliser la console python mais je ne bouge pas :(

Quelqu'un peut-il m'expliquer comment y parvenir pas à pas?

Raquel Ribeiro
la source
voir ma réponse ci-dessous. le problème a déjà été résolu et expliqué
Generic Wevers
2
Et pourquoi ce double avec celui marqué? Peut-être que l'OP essaie d'apprendre les pyqgis et comment utiliser le python si vous considérez ses caractères gras.
nickves
Veuillez préciser votre question. Voulez-vous ne pas les charger manuellement dans QGIS? Voulez-vous les convertir dans un autre format? Quelle est exactement ta question?
bugmenot123
1. Importez tous les fichiers en un seul processus dans qgis 2. projetez-les 3. exportez-les à nouveau en tant que csv mais en coordonnées utm
Raquel Ribeiro
cat * .csv> one_file.csv (ou quel que soit l'équivalent de Windows) combinera tous vos fichiers csv en un seul. 500 n'est vraiment pas un si grand nombre :-)
John Powell

Réponses:

15

Si vous cherchez à reprojeter des fichiers csv à partir de la console Python dans QGIS, vous pouvez utiliser le script suivant. Il vous suffirait de changer les trois chemins mentionnés dans les commentaires.

Essentiellement, le script importe vos fichiers csv dans QGIS en tant que fichiers de formes (en supposant que vos champs géométriques sont nommés Xet Y). Il utilise ensuite les algorithmes qgis:reprojectlayeret qgis:fieldcalculatorde la boîte à outils de traitement pour reprojeter et mettre à jour les champs Xet Yavec les nouvelles coordonnées. Il les enregistre ensuite dans un dossier et les convertit en fichiers csv dans un chemin que vous spécifiez. Donc, à la fin, vous avez mis à jour les fichiers de formes et les fichiers csv dans des dossiers séparés.

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")   

J'espère que cela t'aides!

Joseph
la source
2
Excellente réponse - vous avez tout là-bas!. Une question si cela ne vous dérange pas: vous devez toujours ajouter / supprimer des couches à QgsMapLayerRegistry même si vous faites des choses à partir de la console python?
nickves
1
@nickves - Haha merci beaucoup mon pote! Hmm je pourrais ne pas avoir à ajouter / supprimer des couches (je suis certain que le script peut être considérablement réduit). Je ne suis pas un expert mais je le testerai plus tard et je vous recontacterai. À moins que vous ne puissiez fournir un script beaucoup plus soigné, auquel cas vous devriez le poster comme réponse, je le voterais favorablement :)
Joseph
@nickves - Merci encore pour votre copain de suggestion! Le code a été modifié pour éviter d'ajouter / supprimer les couches une deuxième fois :)
Joseph
@RaquelRibeiro - Bienvenue! Heureux que cela ait été utile :)
Joseph
@Joseph puis-je vous demander quelque chose à nouveau? Aux lignes 14 et 15, les nombres: 0, 10, 10 définissent exactement quoi? (les coordonnées de sortie ont trop de zéros à droite et je veux les minimiser)
Raquel Ribeiro
8

Une solution rapide pour transformer un fichier séparé par espace contenant "lon lat" dans WGS84 en UTM33N mais vous n'obtenez aucune autre donnée:

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

Cela fonctionne et cela préserve l'ordre des données, alors peut-être une autre boucle utilisant par exemple awk pour combiner les données descriptives avec les coordonnées?

Éditer. En raison des commentaires désordonnés que j'ai faits ci-dessous, je modifierai la réponse ici à la place.

Le script suivant devrait faire le travail de lecture de plusieurs fichiers csv, en ajoutant de nouvelles colonnes de coordonnées à chaque fichier.

#!/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

Sur OSX, vous devrez installer la dernière version (2009) de sed et utiliser la première ligne non commentée de la boucle. Pour Linux, commentez le premier et utilisez le second. Ajustez le en -F " "fonction du format du séparateur dans vos fichiers csv, par exemple -F ","pour les virgules séparées. Notez également que la transformation d'élévation se fait vers l'ellipsoïde, pas le géoïde, alors assurez-vous de transformer les hauteurs en conséquence.

mercergeoinfo
la source
Je viens de me souvenir d'avoir fait quelque chose de similaire il y a quelque temps et d'avoir publié une solution sur mon blog. Il est écrit pour Mac mais est basé sur bash. La plus grande différence est le problème avec sed sur OS X, que je traite à la fin du post: mercergeoinfo.blogspot.se/2014/01/…
mercergeoinfo
Le dernier commentaire était un peu désordonné. Utilisez cette ligne dans le script bash ci-dessus pour parcourir tous les fichiers 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}Remplacez / usr / local / sed par juste sed si vous n'êtes pas sous OSX. Ce n'est pas idéal si vos fichiers csv sont séparés par des espaces, comme le suppose la ligne ci-dessus, mais cela fonctionne. Si vous avez séparé la virgule-F " "-F ","
passez
Je me demande pourquoi le code mis à jour est dans les commentaires et vous n'avez pas mis à jour votre réponse ci-dessus. Le code dans le commentaire est vraiment difficile à lire. Voyez-vous le lien d'édition sous votre réponse?
Miro
Oui, mais ce n'était pas vraiment une mise à jour, plutôt une additionnelle. Assez salissant, je suis d'accord. Je suppose que je devrais mettre à jour la réponse d'origine. Merci
mercergeoinfo
7

Utiliser qgis ou même OGR est exagéré pour cela.
Utilisez pyproj( https://pypi.python.org/pypi/pyproj ) combiné avec l'écrivain cython python et quelques astuces de bibliothèque standard. Vous n'avez pas besoin d'installer autre chose que pyprojcela!

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
blord-castillo
la source
Je me rends compte que l'affiche est assez inexpérimentée avec le python. Je n'utilise pas régulièrement QGIS, donc quelqu'un avec plus d'expérience avec cette plateforme pourrait-il expliquer où python est installé? L'affiche devrait en faire un script autonome et probablement l'exécuter depuis IDLE. Je n'ai pas d'installation en cours, donc je ne sais pas si pyprojdoit être installé séparément pour l'affiche, ou est déjà là.
blord-castillo
1
jamais utilisé de fonction partielle auparavant. Fera désormais. +1
nickves
4

Vous n'avez pas besoin de python. Utilisez simplement la ligne de commande et ogr2ogr. Dans votre cas, le plus important est le paramètre -t_srs srs_def.

Ceci est déjà expliqué dans cette réponse à Comment puis-je convertir un fichier Excel avec des colonnes x, y en un fichier de formes?

MISE À JOUR Je n'ai pas le temps de vous écrire votre code complet. Mais le problème sera qu'il a besoin d'un peu plus de code en python que vous ne le pensez.

Votre principal problème sera que travailler avec des fichiers csv n'est pas aussi confortable que d'utiliser des fichiers de formes. Ainsi, vous devrez d'abord convertir le csv en forme qui nécessite un fichier VRT. Ceci est expliqué dans le premier lien. Ici, vous devrez écrire un script python en boucle dans vos fichiers qui génère automatiquement les fichiers vrt.

Ceci est un script que j'ai utilisé moi-même. Vous devez tester si cela fonctionne pour vous. J'ai déjà inclus la conversion de WGS 84 en 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)

Vous devez ajuster les paramètres du nom de champ , src , x et y en fonction de votre fichier csv.

UPDATE2

Après réflexion, je me demande pourquoi voulez-vous utiliser QGIS? Vous pouvez utiliser un script python comme celui-ci pour convertir directement vos coordonnées de WGS en UTM. Dans ce cas, c'est un simple csv ouvert, lire les coordonnées, transformer les coordonnées et l'enregistrer dans un nouveau fichier.

Generic Wevers
la source
Je pense que ce n'est pas ce que je recherche ... J'ai près de 500 fichiers csv (coordonnées wgs84) et voici ce que je veux faire: 1. Importez tous les fichiers csv en une seule fois dans q gis 2. projetez-les 3. les exporter dans des fichiers csv (encore) mais avec des coordonnées différentes (conversion en utm33N)
Raquel Ribeiro
je pense que j'ai besoin d'un processus par lots ou quelque chose comme ça pour le faire ...
Raquel Ribeiro
4
mais pourquoi voulez-vous faire ça? 1. vous pouvez faire la même chose (ce que vous avez décrit) à partir de la ligne de commande sans qgis. 2. vous pouvez le faire en mode batch. 3. en python, c'est presque la même chose. vous utiliseriez également ogr2ogr
Generic Wevers
2
"Simplement" en utilisant la ligne de commande n'est vraiment pas une réponse. La ligne de commande n'est jamais facile à utiliser si vous ne savez pas comment le faire. Et je ne trouve vraiment pas la solution dans la réponse liée. Pourquoi ne pas simplement donner au pauvre garçon un exemple de lot avec ogr2ogr, et tout irait bien?
Bernd V.
1
ok, 1. vous pouvez lire gis.stackexchange.com/help/how-to-ask . après cela et 5 minutes de recherche sur Google, vous admettrez que la question est très mal étudiée et peut être résolue avec des réponses déjà données. 2. Si cela ne peut toujours pas être résolu, je suppose que tout le monde sera heureux de vous aider. mais comme je suis une bonne personne, je vais vous donner quelques conseils supplémentaires.
Generic Wevers