Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

this is probably a very simple question, but I'm very new to GIS so I just need one last nudge. I've got a sense of the Linux utility gdaltransform and I am feeding it coordinate pairs in line-by-line format in "interactive mode" where the program waits for you to input coords and then it just spits back the results on the console.

But a better way for scripting purposes would be able to just feed it an entire polygon of points and have it spit out the re-projected points, all in one shot.

For example, here's what I'd like to do:

gdaltransform -s_srs EPSG:26986 -t_srs EPSG:4326 255771.72999999672,951430.5599999987 255756.3400000036,951443.129999999 255733.53000000122,951443.4400000013 255700.79999999702,951404.4400000015 255693.40999999642,951379.1900000013 255650.22999999672,951370.9400000012 255614.3400000036,951413.879999999 255561.04999999705,951408.3099999985 255510.45000000298,951377.7500000001 255441.9099999964,951361.75 255406.81000000238,951310.6299999988 255409.90999999642,951242.6900000013 255428.39000000063,951161.129999999 255434.14000000063,951070.2500000001 255498.34000000358,950986.0599999987 255507.03000000116,950854.1900000013 255545.59000000358,950798.25 255545.7800000012,950772.129999999 255545.79519999772,950770.4767999984 255558.06970000267,950778.8845000006 255617.48359999806,950816.685800001 255754.21130000052,950904.0976999998 255789.04469999668,950926.2617999986 255798.64760000262,950932.4488000013 255894.92400000244,950993.9884000011 255916.64130000022,951009.1233000011 255916.1400000006,951017.1900000012 255905.78000000116,951188.379999999 255775.32999999818,951303.5599999987 255764.8900000006,951334.3099999987 255771.72999999672,951430.5599999987

Simple, right? Unfortunately, the program reads the beginning of this string as the name of an input file. (So I tried to just dump this list into an input source file and point it to that, but it says the file is in the 'wrong format'.)

How do I get this to work to process bulk pairs, in non-interactive mode?

share|improve this question
add comment

2 Answers

The following should work, but for me the output file is always zero length (on Windows, the output redirection works fine on Linux).

gdaltransform -s_srs EPSG:26986 -t_srs EPSG:4326 < in.txt > out.txt

Where in.txt contain a line separated list of space separated coordinate pairs. If you remove the > out.txt you'll see the list of transformed coordinate pairs printed to the terminal.

The file in 'wrong format' error message is because gdaltransform is expecting projection definition files when filenames are given on the command line.

Edit: below is a python script that will reproject a csv file:

from osgeo import osr,ogr
import sys,csv

class Transform(object):
    def __init__(self,s_srs,t_srs,in_csv,out_csv,delimiter=',',quotechar='|'):
        self.s_srs = osr.SpatialReference()
        self.t_srs = osr.SpatialReference()
        self.s_srs.ImportFromEPSG(s_srs)
        self.t_srs.ImportFromEPSG(t_srs)
        infile = open(in_csv, 'rb')
        outfile = open(out_csv, 'wb')
        csvreader = csv.reader(infile, delimiter=delimiter, quotechar=quotechar)
        csvwriter = csv.writer(outfile, delimiter=delimiter, quotechar=quotechar)
        for row in csvreader:
            inx,iny=[float(coord) for coord in row]
            csvwriter.writerow(self.transform(inx,iny))

        del csvreader
        del csvwriter
        infile.close()
        outfile.close()

    def transform(self,xcoord,ycoord):
        geom = ogr.Geometry(ogr.wkbPoint)
        geom.SetPoint_2D(0, xcoord,ycoord)
        geom.AssignSpatialReference(self.s_srs)
        geom.TransformTo(self.t_srs)
        return geom.GetPoint_2D()

if __name__=='__main__':
    #testing
    #s_srs=26986
    #t_srs=4326
    #in_csv='in_file.csv'
    #out_csv='out_file.csv'

    s_srs=int(sys.argv[1])
    t_srs=int(sys.argv[2])
    in_csv=sys.argv[3]
    out_csv=sys.argv[4]

    Transform(s_srs,t_srs,in_csv,out_csv)
share|improve this answer
    
You shouldn't need the redirection for the input, by default gdaltransform takes an input file. –  kyle Oct 10 '12 at 15:04
    
It seems to only take a gdal supported dataset as input.Line 209 of gdaltransform.cpp: hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly ); –  Luke Oct 10 '12 at 20:06
    
Ah, yes. I was thinking of cs2cs, which gdaltransform wraps. My mistake. –  kyle Oct 10 '12 at 20:50
add comment

to convert a txt file with a list of coordinates from one CRS to another the right tools is cs2cs utility: http://proj.maptools.org/man_cs2cs.html

You find it also in gdal distribution.

The syntax is something like:

cs2cs +init=epsg:32633 +to +init=epsg:4326 c:\tmp\input.txt > c:\tmp\output.txt

It's not probably what you want, but it should be useful for you.

Cheers,

Andrea

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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