3 # Glyn Astill - 11/10/2012
4 # Script to automate cs2cs transform from esatings/northings to latitudes/longitudes
10 use Getopt::Long qw/GetOptions/;
12 use constant false => 0;
13 use constant true => 1;
15 my $csv = Text::CSV->new();
17 my $usage = '-i <codepoint csv file input> -o <latlong csv file output>';
24 my $start_run = time();
27 my $pipe_cs2cs = true;
31 # I'm not positive that the below is setting a the Bursa Wolf or another parameter correctly, hence loss in accuracy.
32 # it's possible to set the parameters individually i.e. "Geo::Proj4->new(proj => "tmerc", ellps => "airy", lat_0 => -49)" which may solve
33 $proj_ng = Geo::Proj4->new("-f '%.7f' +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +no_defs +to +proj=latlong +ellps=WGS84 +towgs84=0,0,0 +nodefs")
36 # Handle command line options
37 Getopt::Long::Configure('no_ignore_case');
39 die $usage unless GetOptions(\%opt, 'infile|i=s', 'outfile|o=s', ) and keys %opt and ! @ARGV;
41 if (!defined($opt{infile})) {
42 print("Please specify an input file.\n");
46 $infile = $opt{infile};
48 if (!defined($opt{outfile})) {
49 print("Please specify an output file.\n");
53 $outfile = $opt{outfile};
56 # Read the codepoint file with eastings and northings, convert and output to outfile
57 # Open our output file which will contain the eastings/northings and longitudes/latitudes and write a headder
58 if (open(OUTCSV, ">", $outfile)) {
59 print OUTCSV ("postcode,easting,northing,latitude,longitude,country_code,admin_county_code,admin_district_code,admin_ward_code\n");
61 # get a list of our files
62 @infiles = glob($infile);
64 foreach my $currfile (@infiles) {
67 print ("Processing $currfile ..");
68 unless (open(INCSV, "<", $currfile)) {
69 print("ERROR: Could not open file:" . $! . ".\n");
72 if ($csv->parse($_)) {
73 @eastnorth=$csv->fields();
77 open(CS2CS,"echo $eastnorth[2] $eastnorth[3] | cs2cs -f '%.7f' +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 +units=m +no_defs +to +proj=latlong +ellps=WGS84 +towgs84=0,0,0 +nodefs |");
80 @latlong = $proj_ng->inverse($eastnorth[2], $eastnorth[3]);
84 @latlong = split(' ', <CS2CS>);
86 print OUTCSV ("$eastnorth[0],$eastnorth[2],$eastnorth[3],$latlong[1],$latlong[0],$eastnorth[4],$eastnorth[7],$eastnorth[8],$eastnorth[9]\n");
89 print OUTCSV ("$eastnorth[0],$eastnorth[2],$eastnorth[3],$latlong[0],$latlong[1],$eastnorth[4],$eastnorth[7],$eastnorth[8],$eastnorth[9]\n");
96 print("ERROR: Could not run command:" . $! . " Line = " . $_ . "\n");
97 print("$eastnorth[2] $eastnorth[3]");
101 if (($linecount%1000) == 0) {
102 print ("..$linecount");
106 $error = $csv->error_input;
107 print("\tFailed to parse line: $error\n");
117 print("ERROR: Could not open file:" . $! . ".\n");
121 $run_time = (($end_run-$start_run)/60);
123 print "Conversion took $run_time minutes\n";
125 print "To run quicker use Geo:proj4 (appears to be less accurate - a bug in my usage?) or stream the data directly to cs2cs in one go rather than constantly calling and opening the commands output\n";