11

I know, my question is similar to some old ones on this site.

I've a lot of CSV files (geo coordinates) to import to qgis (and then to convert them), and the usual way is not the best way to do it (too long).

I have almost 500 CSV files (wgs84 coordinates) and this is what I want to do:

  1. Import all CSV files at once into QGIS
  2. Project them
  3. Export them into CSV files (again) but with different coordinates (convertion to UTM33N)

I'm trying to understand how to use the python console but I'm not moving on :(

Can anyone explain to me how to achieve it step by step?

  • see my answer below. the problem was already solved and explained – Generic Wevers Oct 28 '15 at 11:06
  • 2
    And why is that duplicate with the marked one? Maybe the OP tries to learn pyqgis and how to use python if you consider his/her bolds. – nickves Oct 28 '15 at 11:41
  • Please specify your question. Do you want to not load them manually into QGIS? Do you want to convert them into another format? What exactly is your question? – bugmenot123 Oct 28 '15 at 14:06
  • 1. Import all the files in one process to qgis 2. project them 3. export all them again as csv but in utm coordinates – Raquel Ribeiro Oct 29 '15 at 0:59
  • cat *.csv > one_file.csv (or whatever the windows equivalent is) will combine all your csv files to one. 500 really isn't such a big number :-) – John Powell Oct 29 '15 at 9:07
15

If you're looking to reproject csv files from the Python Console in QGIS then you could use the following script. All you would need to change are the three paths which are mentioned in the comments.

Essentially, the script imports your csv files into QGIS as shapefiles (assuming your geometric fields are named X and Y). It then uses the qgis:reprojectlayer and qgis:fieldcalculator algorithms from the Processing Toolbox to reproject and update the X and Y fields with the new coordinates. It then saves these in a folder and converts them to csv files in a path you specify. So in the end, you have updated shapefiles and csv files in separate folders.

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

Hope this helps!

  • 2
    Great answer - you have it all there!. One question if you don't mind: You still have to add/remove layers at QgsMapLayerRegistry even if you do things from the python console? – nickves Nov 1 '15 at 11:22
  • 1
    @nickves - Haha many thanks buddy! Hmm I might not have to add/remove layers (I'm certain the script can be reduced dramatically). I'm no expert but I will test it later and get back to you. Unless you can provide a much neater script in which case you should post it as an answer, I would upvote it :) – Joseph Nov 2 '15 at 11:01
  • @nickves - Thanks again for your suggestion buddy! Code has been edited to avoid adding/removing the layers a second time :) – Joseph Nov 3 '15 at 11:37
  • @RaquelRibeiro - Most welcome! Glad it was helpful :) – Joseph Nov 10 '15 at 10:43
  • @Joseph may i ask u something again? At lines 14&15, the numbers: 0, 10, 10 are defining exactly what? (the output coordinates have too much zeros on right and i want to minimize them) – Raquel Ribeiro Nov 11 '15 at 1:02
8

A quick solution for transforming a space separated file containing "lon lat" in WGS84 to UTM33N but you don't get any other data:

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

That works and it preserves the order of the data so maybe another loop using e.g. awk to combine the descriptive data with the coordinates?

Edit. Due to the messy comments I made below I'll edit the answer here instead.

The following script should do the job of reading multiple csv files, adding new coordinate columns to each file.

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

On OSX you will need to install the latest (2009) version of sed and use the first, uncommented line in the loop. For Linux comment out the first and use the second. Adjust the -F " " according to the format of the separator in your csv files e.g. -F "," for comma separated. Also note that the elevation transformation is to the ellipsoid,not the geoid, so be sure to transform the heights accordingly.

  • I just remembered doing something similar a while ago and posting a solution to my blog. It's written for Mac but is bash based. The biggest difference is the issue with sed on OS X, which I deal with at the end of the post: mercergeoinfo.blogspot.se/2014/01/… – mercergeoinfo Nov 2 '15 at 9:12
  • The last comment was a bit messy. Use this line in the above bash script to loop through all files 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} Replace the /usr/local/sed with just sed if you aren't on OSX. This isn't ideal if your csv files are space separated, as the above line assumes, but it works. If you have comma separated then change -F " " to -F "," – mercergeoinfo Nov 2 '15 at 9:50
  • I wonder why is the updated code in comments and you did not updated your answer above. The code in the comment is really hard to read. Do you see the edit link below your answer? – Miro Nov 3 '15 at 23:21
  • Yup, but then it wasn't really an update, more like an additional. Pretty messy, I agree. I guess I should update the original answer. Thanks – mercergeoinfo Nov 4 '15 at 8:42
7

Using qgis or even OGR is overkill for this.
Use pyproj (https://pypi.python.org/pypi/pyproj) combined with the python csv writer and a few standard library tricks. You do not need to install anything other than pyproj for this!

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
  • I realize the poster is pretty inexperienced with python. I do not regularly use QGIS, so could someone with more experience with that platform explain where python is installed? The poster should make this a standalone script and probably run it from IDLE. I do not have a current installation, so I don't know if pyproj needs to be installed separately for the poster, or is already there. – blord-castillo Oct 30 '15 at 16:14
  • 1
    never used partial function before. Will do from now on. +1 – nickves Nov 1 '15 at 11:26
4

You don't need python. Simply use the command line and ogr2ogr. In your case most important is the -t_srs srs_def parameter.

This is already explained in this answer to How can I convert an excel file with x, y columns to a shapefile?

UPDATE I don't have the time to write you your complete code. But the problem will be that it needs a little more code in python than you may think.

Your main problem will be that working with csv files is not as comfortable as using shapefiles. Thus you will first need to convert the csv to shape which needs VRT file. This is explained in the first link. Here you will need to write a python script looping through your files which automatically generates the vrt files.

This is a script I used myself. You have to test if it works for you. I already included the conversion from WGS 84 to 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)

You need to adjust the parameters for Field name, src, x and y according to your csv file.

UPDATE2

After some thinking, I ask myself why do you want to use QGIS at all? You could use a python script like this to directly convert your coordinates from WGS to UTM. In this case it is a simple open csv, read coordinate, transform coordinate and save it to a new file.

  • i think this is not what i'm looking for... I have almost 500 csv files (wgs84 coordinates) and this is what i want to do: 1. Import all csv files in once to q gis 2. project them 3. export them into csv files (again) but with diferent coordinates (convertion to utm33N) – Raquel Ribeiro Oct 28 '15 at 11:38
  • i think i need a batch process or somethig like that to do it... – Raquel Ribeiro Oct 28 '15 at 11:46
  • 4
    but why do you want to do that? 1. you can do the same (what you described) from the command line without qgis. 2. you can do this in batch mode. 3. in python it is nearly the same. you would also use ogr2ogr – Generic Wevers Oct 28 '15 at 11:47
  • 2
    "Simply" using the command line is really not an answer. Command line is never easy to use if you have no idea how to do it. And I really can't find the solution in the linked answer. Why not just giving the poor fellow an example batch with ogr2ogr, and everything would be fine? – Bernd V. Oct 29 '15 at 0:00
  • 1
    ok, 1. you may read gis.stackexchange.com/help/how-to-ask. after that and 5 minutes of googleing you will admit, that the question is very poorly researched and can be solved with already given answers. 2. If it still can not be solved, i guess everybody will be pleased to help. but as i am a good person, i will give some more hints. – Generic Wevers Oct 29 '15 at 8:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.