Como reprojetar 500 arquivos CSV de forma eficiente e fácil usando o QGIS?

11

Eu sei, minha pergunta é semelhante a algumas antigas deste site.

Eu tenho muitos arquivos CSV (coordenadas geográficas) para importar para o qgis (e depois convertê-los), e a maneira usual não é a melhor maneira de fazê-lo (muito tempo).

Eu tenho quase 500 arquivos CSV (coordenadas wgs84) e é isso que eu quero fazer:

  1. Importe todos os arquivos CSV de uma vez para o QGIS
  2. Projete-os
  3. Exporte-os para arquivos CSV (novamente), mas com coordenadas diferentes (conversão para UTM33N)

Estou tentando entender como usar o console python, mas não estou seguindo em frente :(

Alguém pode me explicar como alcançá-lo passo a passo?

Raquel Ribeiro
fonte
veja minha resposta abaixo. o problema já foi resolvido e explicado
Generic Wevers
2
E por que isso é duplicado com o marcado? Talvez o OP tente aprender pyqgis e como usá-lo se você considerar seus negritos.
Nick28 #
Especifique sua pergunta. Deseja não carregá-los manualmente no QGIS? Deseja convertê-los em outro formato? Qual é a tua pergunta exatamente?
precisa
1. Importe todos os arquivos em um único processo para qgis 2. projeto los 3. exportar todos eles novamente como CSV, mas em coordenadas UTM
Raquel Ribeiro
cat * .csv> one_file.csv (ou qualquer que seja o equivalente do Windows) combinará todos os seus arquivos csv em um. 500 realmente não é um número tão grande :-)
John Powell

Respostas:

15

Se você deseja reprojetar arquivos csv do Python Console no QGIS, use o seguinte script. Tudo que você precisa mudar são os três caminhos mencionados nos comentários.

Essencialmente, o script importa seus arquivos csv no QGIS como shapefiles (assumindo que seus campos geométricos sejam nomeados Xe Y). Em seguida, ele usa os algoritmos qgis:reprojectlayere qgis:fieldcalculatorda Caixa de Ferramentas de Processamento para reprojetar e atualizar os campos Xe Ycom as novas coordenadas. Em seguida, os salva em uma pasta e os converte em arquivos csv no caminho especificado. Portanto, no final, você atualizou os arquivos shapefiles e os arquivos csv em pastas 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 isto ajude!

Joseph
fonte
2
Ótima resposta - você tem tudo lá! Uma pergunta se você não se importa: você ainda precisa adicionar / remover camadas no QgsMapLayerRegistry, mesmo que faça coisas no console python?
nickves
1
@nickves - Haha muito obrigado amigo! Hmm, talvez não seja necessário adicionar / remover camadas (tenho certeza de que o script pode ser reduzido drasticamente). Não sou especialista, mas testarei mais tarde e retornarei a você. A menos que você pode fornecer um roteiro mais puro tanto no caso em que você deve publicá-la como uma resposta, gostaria de upvote-lo :)
Joseph
@ Nickves - Mais uma vez obrigado pelo seu amigo sugestão! Código foi editado para evitar adicionar / remover as camadas uma segunda vez :)
Joseph
@RaquelRibeiro - Muito bem-vindo! Ainda bem que foi :) útil
Joseph
@ Joseph posso perguntar-lhe algo de novo? Nas linhas 14 e 15, os números: 0, 10, 10 estão definindo exatamente o que? (as coordenadas de saída têm muito zeros à direita e eu quero minimizá-los)
Raquel Ribeiro
8

Uma solução rápida para transformar um arquivo separado por espaço que contém "long lat" em WGS84 em UTM33N, mas você não obtém outros dados:

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

Isso funciona e preserva a ordem dos dados, então talvez outro loop usando, por exemplo, awk para combinar os dados descritivos com as coordenadas?

Editar. Devido aos comentários confusos que fiz abaixo, editarei a resposta aqui.

O script a seguir deve fazer o trabalho de ler vários arquivos csv, adicionando novas colunas de coordenadas a cada arquivo.

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

No OSX, você precisará instalar a versão mais recente do sed (2009) e usar a primeira linha não comentada do loop. Para o Linux, comente o primeiro e use o segundo. Ajuste de -F " "acordo com o formato do separador nos arquivos csv, por exemplo, -F ","para separar vírgulas. Observe também que a transformação de elevação é no elipsóide, não no geóide, portanto, certifique-se de transformar as alturas de acordo.

mercergeoinfo
fonte
Acabei de me lembrar de fazer algo semelhante há um tempo atrás e postar uma solução no meu blog. Ele foi escrito para Mac, mas é baseado em bash. A maior diferença é o problema com o sed no OS X, com o qual lido no final do post: mercergeoinfo.blogspot.se/2014/01/…
mercergeoinfo
O último comentário foi um pouco confuso. Use esta linha no script bash acima para percorrer todos os arquivos 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}Substitua / usr / local / sed por just sed se você não estiver no OSX. Isso não é ideal se os arquivos csv estiverem separados por espaço, como supõe a linha acima, mas funciona. Se você tiver separados por vírgula, mude -F " "para-F ","
mercergeoinfo
Gostaria de saber por que o código atualizado está nos comentários e você não atualizou sua resposta acima. O código no comentário é realmente difícil de ler. Você vê o link de edição abaixo da sua resposta?
Miro
Sim, mas não foi realmente uma atualização, mais como uma adicional. Muito bagunçado, eu concordo. Acho que devo atualizar a resposta original. Obrigado
mercergeoinfo 4/15
7

Usar qgis ou mesmo OGR é um exagero para isso.
Use pyproj( https://pypi.python.org/pypi/pyproj ) combinado com o gravador python csv e alguns truques padrão da biblioteca. Você não precisa instalar nada além pyprojdisso!

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
fonte
Sei que o pôster é inexperiente com python. Como não uso o QGIS regularmente, alguém poderia ter mais experiência com essa plataforma para explicar onde o python está instalado? O pôster deve torná-lo um script independente e provavelmente executá-lo a partir do IDLE. Como não tenho uma instalação atual, não sei se pyprojprecisa ser instalado separadamente para o pôster ou se já está lá.
Blord-castillo 30/10/2015
1
nunca usei função parcial antes. Vai fazer a partir de agora. +1
nickves 01/11/2015
4

Você não precisa de python. Basta usar a linha de comando e ogr2ogr. No seu caso, o mais importante é o parâmetro -t_srs srs_def.

Isso já está explicado nesta resposta para Como converter um arquivo do Excel com colunas x, y em um shapefile?

ATUALIZAÇÃO Não tenho tempo para escrever seu código completo. Mas o problema será que ele precisa de um pouco mais de código em python do que você imagina.

Seu principal problema será que trabalhar com arquivos csv não é tão confortável quanto usar shapefiles. Assim, você primeiro precisará converter o csv para a forma que precisa do arquivo VRT. Isso é explicado no primeiro link. Aqui, você precisará escrever um script python em loop pelos arquivos, que gera automaticamente os arquivos vrt.

Este é um script que eu me usei. Você tem que testar se funciona para você. Eu já incluí a conversão de WGS 84 para 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)

Você precisa ajustar os parâmetros para Nome do campo , src , x e y de acordo com o seu arquivo csv.

UPDATE2

Depois de pensar um pouco, pergunto-me por que você deseja usar o QGIS? Você pode usar um script python como este para converter diretamente suas coordenadas de WGS em UTM. Nesse caso, é um csv aberto simples, lê coordenadas, transforma coordenadas e salva-as em um novo arquivo.

Wevers genéricos
fonte
Eu acho que não é isso que estou procurando ... Tenho quase 500 arquivos csv (coordenadas wgs84) e é isso que eu quero fazer: 1. Importe todos os arquivos csv de uma vez para q gis 2. projete-os 3. exportá-los para arquivos CSV (de novo), mas com coordenadas diferentes (conversão para utm33N)
Raquel Ribeiro
eu acho que eu preciso de um processo em lote ou somethig assim fazê-lo ...
Raquel Ribeiro
4
mas por que você quer fazer isso? 1. você pode fazer o mesmo (o que você descreveu) na linha de comando sem o qgis. 2. você pode fazer isso no modo de lote. 3. em python, é quase o mesmo. você também usaria ogr2ogr
Generic Wevers
2
"Simplesmente" usar a linha de comando não é realmente uma resposta. A linha de comando nunca é fácil de usar se você não tem idéia de como fazê-lo. E realmente não consigo encontrar a solução na resposta vinculada. Por que não apenas dar ao pobre sujeito um lote de exemplo com ogr2ogr, e tudo ficaria bem?
Bernd V.
1
ok, 1. você pode ler gis.stackexchange.com/help/how-to-ask . depois disso e 5 minutos de pesquisa no Google, você admitirá que a pergunta é muito pouco pesquisada e pode ser resolvida com respostas já fornecidas. 2. Se ainda não puder ser resolvido, acho que todos terão prazer em ajudar. mas como sou uma boa pessoa, darei mais algumas dicas.
Wevers genéricos