Как эффективно и легко перепроектировать 500 файлов CSV с помощью QGIS?

11

Я знаю, мой вопрос похож на некоторые старые на этом сайте.

У меня много CSV-файлов (географических координат) для импорта в qgis (а затем для их преобразования), и обычный способ не самый лучший способ сделать это (слишком долго).

У меня есть почти 500 файлов CSV (координаты wgs84), и это то, что я хочу сделать:

  1. Импортируйте все файлы CSV одновременно в QGIS
  2. Проецируйте их
  3. Экспортируйте их в файлы CSV (снова), но с другими координатами (преобразование в UTM33N)

Я пытаюсь понять, как использовать консоль Python, но я не двигаюсь дальше :(

Может кто-нибудь объяснить мне, как это сделать шаг за шагом?

Ракель Рибейро
источник
см. мой ответ ниже. проблема была уже решена и объяснена
Generic Wevers
2
И почему этот дубликат с отмеченным? Может быть, ОП пытается выучить Pyqgis и как использовать Python, если вы считаете его / ее смелыми.
Никв
Пожалуйста, уточните свой вопрос. Вы не хотите загружать их вручную в QGIS? Вы хотите конвертировать их в другой формат? Что именно ваш вопрос?
bugmenot123
1. Импортируйте все файлы в одном процессе в qgis 2. Спроецируйте их 3. Экспортируйте все снова как csv, но в координатах utm
Ракель Рибейро
cat * .csv> one_file.csv (или любой другой эквивалент Windows) объединит все ваши CSV-файлы в один. 500 действительно не такое большое число :-)
Джон Пауэлл

Ответы:

15

Если вы хотите перепроектировать CSV-файлы из консоли Python в QGIS, вы можете использовать следующий скрипт. Все, что вам нужно изменить, это три пути, которые упоминаются в комментариях.

По сути, скрипт импортирует ваши CSV-файлы в QGIS как шейп-файлы (при условии, что ваши геометрические поля названы Xи Y). Затем он использует qgis:reprojectlayerи qgis:fieldcalculatorалгоритмы из Processing Toolbox для перепроецировать и обновления Xи Yполя с новыми координатами. Затем он сохраняет их в папке и преобразует их в CSV-файлы по указанному вами пути. Итак, в конце вы обновили шейп-файлы и CSV-файлы в отдельных папках.

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

Надеюсь это поможет!

Джозеф
источник
2
Отличный ответ - у вас все это есть! Один вопрос, если вы не возражаете: вам все еще нужно добавлять / удалять слои в QgsMapLayerRegistry, даже если вы делаете что-то с консоли Python?
nickves
1
@nickves - Хаха, большое спасибо, приятель! Хм, возможно, мне не придется добавлять / удалять слои (я уверен, что сценарий может быть значительно сокращен). Я не эксперт, но позже опробую и свяжусь с вами. Если вы не можете предоставить более аккуратный сценарий, в этом случае вы должны опубликовать его в качестве ответа, я бы высказался за него :)
Joseph
@nickves - Еще раз спасибо за ваше предложение приятель! Код был отредактирован, чтобы избежать добавления / удаления слоев во второй раз :)
Joseph
@RaquelRibeiro - Добро пожаловать! Рад, что это было полезно :)
Иосиф
@ Джозеф, могу я еще кое-что спросить? В строках 14 и 15 числа: 0, 10, 10 определяют, что именно? (выходные координаты имеют слишком много нулей справа, и я хочу минимизировать их)
Ракель Рибейро
8

Быстрое решение для преобразования файла, разделенного пробелом, содержащего «lon lat» в WGS84, в UTM33N, но вы не получите никаких других данных:

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

Это работает, и это сохраняет порядок данных, так что, может быть, другой цикл, используя, например, awk, чтобы объединить описательные данные с координатами?

Редактировать. Из-за грязных комментариев, которые я сделал ниже, я вместо этого отредактирую ответ.

Следующий скрипт должен выполнять чтение нескольких файлов CSV, добавляя новые столбцы координат в каждый файл.

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

В OSX вам нужно будет установить последнюю версию sed (2009) и использовать первую некомментированную строку в цикле. Для Linux закомментируйте первое и используйте второе. Отрегулируйте в -F " "соответствии с форматом разделителя в ваших CSV-файлах, например, через -F ","запятую. Также обратите внимание, что преобразование высот относится к эллипсоиду, а не к геоиду, поэтому обязательно измените высоту соответствующим образом.

mercergeoinfo
источник
Я только что вспомнил, как делал нечто подобное некоторое время назад и публиковал решение для своего блога. Он написан для Mac, но основан на bash. Самым большим отличием является проблема с sed на OS X, с которой я сталкиваюсь в конце поста: mercergeoinfo.blogspot.se/2014/01/…
mercergeoinfo
Последний комментарий был немного грязным. Используйте эту строку в приведенном выше скрипте bash для циклического перебора всех файлов. 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}Замените / usr / local / sed на просто sed, если вы не используете OSX. Это не идеально, если ваши CSV-файлы разделены пробелом, как предполагает приведенная выше строка, но это работает. Если вы разделяете запятую, измените -F " "на-F ","
mercergeoinfo
Интересно, почему обновленный код в комментариях, а вы не обновили свой ответ выше. Код в комментарии действительно трудно читать. Вы видите ссылку для редактирования ниже вашего ответа?
Миро
Да, но тогда это было не совсем обновление, больше похоже на дополнительное. Довольно грязно, я согласен. Я думаю, что я должен обновить оригинальный ответ. Спасибо
mercergeoinfo
7

Использование qgis или даже OGR для этого излишне.
Используйте pyproj( https://pypi.python.org/pypi/pyproj ) в сочетании с Python Csv Writer и несколько стандартных библиотечных приемов. Вам не нужно устанавливать ничего, кроме как pyprojдля этого!

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
источник
Я понимаю, что плакат довольно неопытен с питоном. Я не пользуюсь QGIS регулярно, поэтому кто-нибудь с большим опытом работы с этой платформой может объяснить, где установлен python? Постер должен сделать это автономным скриптом и, вероятно, запустить его из IDLE. У меня нет текущей установки, поэтому я не знаю, pyprojнужно ли устанавливать отдельно для плаката или уже есть.
blord-castillo
1
никогда раньше не использовал частичную функцию. Будет делать с этого момента. +1
ники
4

Вам не нужен питон. Просто используйте командную строку и ogr2ogr. В вашем случае наиболее важным является параметр -t_srs srs_def.

Это уже объясняется в этом ответе на вопрос Как я могу преобразовать файл Excel с столбцами x, y в шейп-файл?

ОБНОВЛЕНИЕ У меня нет времени, чтобы написать вам ваш полный код. Но проблема будет в том, что ему нужно немного больше кода на python, чем вы думаете.

Ваша главная проблема будет в том, что работать с CSV-файлами не так удобно, как с использованием шейп-файлов. Таким образом, вам сначала нужно будет преобразовать CSV в форму, для которой нужен файл VRT. Это объясняется в первой ссылке. Здесь вам нужно написать скрипт Python, который просматривает ваши файлы и автоматически генерирует файлы vrt.

Это сценарий, который я использовал сам. Вы должны проверить, работает ли это для вас. Я уже включил преобразование из WGS 84 в 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)

Вам необходимо настроить параметры для поля Имя , src , x и y в соответствии с вашим CSV-файлом.

UPDATE2

Подумав немного, я спрашиваю себя, почему вы вообще хотите использовать QGIS? Вы можете использовать питон скрипт , как это напрямую конвертировать ваши координаты из WGS в UTM. В этом случае это просто открыть CSV, прочитать координаты, преобразовать координаты и сохранить его в новый файл.

Общие Wevers
источник
я думаю, что это не то, что я ищу ... У меня есть почти 500 CSV-файлов (координаты WGS84), и это то, что я хочу сделать: 1. Импортировать все CSV-файлы за один раз в Q-GIS 2. Проектировать их 3. экспортируйте их в CSV-файлы (снова), но с разными координатами (преобразование в utm33N)
Ракель Рибейро
я думаю, что мне нужен пакетный процесс или что-то подобное, чтобы сделать это ...
Ракель Рибейро
4
но почему ты хочешь это сделать? 1. вы можете сделать то же самое (что вы описали) из командной строки без qgis. 2. Вы можете сделать это в пакетном режиме. 3. в питоне это почти то же самое. Вы также использовали бы ogr2ogr
Generic Wevers
2
«Просто» с помощью командной строки на самом деле не является ответом. Командная строка никогда не бывает простой в использовании, если вы не знаете, как это сделать. И я действительно не могу найти решение в связанном ответе. Почему бы просто не дать беднягу пример пакета с ogr2ogr, и все будет хорошо?
Бернд В.
1
хорошо, 1. вы можете прочитать gis.stackexchange.com/help/how-to-ask . после этого и 5 минут гугления вы признаете, что вопрос очень плохо исследован и может быть решен с помощью уже предоставленных ответов. 2. Если это все еще не может быть решено, я думаю, что все будут рады помочь. но так как я хороший человек, я дам еще несколько советов.
Общие Wevers