git.8kb.co.uk Git - postgresql/geographic_data/blob - os_convert.pl

   1 #!/usr/bin/perl
   2
   3 # Glyn Astill - 11/10/2012
   4 # Script to automate cs2cs transform from esatings/northings to latitudes/longitudes
   5
   6 use strict;
   7 use warnings;
   8 use Text::CSV;
   9 use Geo::Proj4;
  10 use Getopt::Long qw/GetOptions/;
  11
  12 use constant false => 0;
  13 use constant true  => 1;
  14
  15 my $csv = Text::CSV->new();
  16 my $error;
  17 my $usage = '-i <codepoint csv file input> -o <latlong csv file output>';
  18 my $infile;
  19 my $outfile;
  20 my @latlong;
  21 my @eastnorth;
  22 my $linecount;
  23 my @infiles;
  24 my $start_run = time();
  25 my $end_run;
  26 my $run_time;
  27 my $pipe_cs2cs = true;
  28 my $proj_ng;
  29
  30 if (!$pipe_cs2cs) {
  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")
  34 }
  35
  36 # Handle command line options
  37 Getopt::Long::Configure('no_ignore_case');
  38 use vars qw{%opt};
  39 die $usage unless GetOptions(\%opt, 'infile|i=s', 'outfile|o=s', ) and keys %opt and ! @ARGV;
  40
  41 if (!defined($opt{infile})) {
  42         print("Please specify an input file.\n");
  43         die $usage;
  44 }
  45 else {
  46         $infile = $opt{infile};
  47 }
  48 if (!defined($opt{outfile})) {
  49         print("Please specify an output file.\n");
  50         die $usage;
  51 }
  52 else {
  53         $outfile = $opt{outfile};
  54 }
  55
  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");
  60
  61         # get a list of our files
  62         @infiles = glob($infile);
  63
  64         foreach my $currfile (@infiles) {
  65                 # read in the file
  66                 $linecount = 0;
  67                 print ("Processing $currfile ..");
  68                 unless (open(INCSV, "<", $currfile)) {
  69                         print("ERROR: Could not open file:" . $! . ".\n");
  70                 }
  71                 while (<INCSV>) {
  72                         if ($csv->parse($_)) {
  73                                 @eastnorth=$csv->fields();
  74
  75                                 eval {
  76                                         if ($pipe_cs2cs) {
  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 |");
  78                                         }
  79                                         else {
  80                                                 @latlong = $proj_ng->inverse($eastnorth[2], $eastnorth[3]);
  81                                         }
  82
  83                                         if ($pipe_cs2cs) {
  84                                                 @latlong =  split(' ', <CS2CS>);
  85                                                 close CS2CS;
  86                                                 print OUTCSV ("$eastnorth[0],$eastnorth[2],$eastnorth[3],$latlong[1],$latlong[0],$eastnorth[4],$eastnorth[7],$eastnorth[8],$eastnorth[9]\n");
  87                                         }
  88                                         else {
  89                                                 print OUTCSV ("$eastnorth[0],$eastnorth[2],$eastnorth[3],$latlong[0],$latlong[1],$eastnorth[4],$eastnorth[7],$eastnorth[8],$eastnorth[9]\n");
  90
  91
  92
  93                                         }
  94                                 };
  95                                 if ($@) {
  96                                         print("ERROR: Could not run command:" . $! . " Line = " . $_ . "\n");
  97                                         print("$eastnorth[2] $eastnorth[3]");
  98                                 }
  99
 100                                 $linecount++;
 101                                 if (($linecount%1000) == 0) {
 102                                         print ("..$linecount");
 103                                 }
 104                         }
 105                         else {
 106                                 $error = $csv->error_input;
 107                                 print("\tFailed to parse line: $error\n");
 108                         }
 109                 }
 110                 close (INCSV);
 111                 print ("..OK\n");
 112         }
 113
 114         close (OUTCSV);
 115 }
 116 else {
 117         print("ERROR: Could not open file:" . $! . ".\n");
 118 }
 119
 120 $end_run = time();
 121 $run_time = (($end_run-$start_run)/60);
 122
 123 print "Conversion took $run_time minutes\n";
 124 if ($pipe_cs2cs) {
 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";
 126 }
 127
 128 exit(0);