Add postcode importer and converter scripts

scripts/OSGB36.rb

@@ -0,0 +1,181 @@
+# https://raw.github.com/LambethCouncil/OSGB36_Converter/master/OSGB36.rb
+module OSGB36
+
+  extend self
+
+  # Takes OSGB36 Easting/Northing coords 
+  # and returns WGS84 Latitude and Longitude
+  def en_to_ll(easting, northing)
+    @OSGB36_ll = to_OSGB36(easting, northing)
+    to_WGS84(@OSGB36_ll[:latitude], @OSGB36_ll[:longitude])
+  end
+
+  # Takes OSGB36 Easting/Northing coords and returns
+  # OSGB36 Latitude and Longitude
+  def to_OSGB36(easting, northing)
+  
+       @OSGB_F0  = 0.9996012717
+       @N0       = -100000.0
+       @E0       = 400000.0
+       @phi0     = deg_to_rad(49.0)
+       @lambda0  = deg_to_rad(-2.0)
+       @a        = 6377563.396
+       @b        = 6356256.909
+       @eSquared = ((@a * @a) - (@b * @b)) / (@a * @a)
+       @phi      = 0.0
+       @lambda   = 0.0
+       @E        = easting
+       @N        = northing
+       @n        = (@a - @b) / (@a + @b)
+       @M        = 0.0
+       @phiPrime = ((@N - @N0) / (@a * @OSGB_F0)) + @phi0
+       
+       begin
+         
+         @M =
+           (@b * @OSGB_F0)\
+             * (((1 + @n + ((5.0 / 4.0) * @n * @n) + ((5.0 / 4.0) * @n * @n * @n))\
+               * (@phiPrime - @phi0))\
+               - (((3 * @n) + (3 * @n * @n) + ((21.0 / 8.0) * @n * @n * @n))\
+                 * Math.sin(@phiPrime - @phi0)\
+                 * Math.cos(@phiPrime + @phi0))\
+               + ((((15.0 / 8.0) * @n * @n) + ((15.0 / 8.0) * @n * @n * @n))\
+                 * Math.sin(2.0 * (@phiPrime - @phi0))\
+                 * Math.cos(2.0 * (@phiPrime + @phi0)))\
+               - (((35.0 / 24.0) * @n * @n * @n)\
+                 * Math.sin(3.0 * (@phiPrime - @phi0))\
+                 * Math.cos(3.0 * (@phiPrime + @phi0))))
+                 
+         @phiPrime += (@N - @N0 - @M) / (@a * @OSGB_F0)
+       
+       end while ((@N - @N0 - @M) >= 0.001)
+       
+       @v = @a * @OSGB_F0 * ((1.0 - @eSquared * sin_pow_2(@phiPrime)) ** -0.5)
+       @rho =
+         @a\
+           * @OSGB_F0\
+           * (1.0 - @eSquared)\
+           * ((1.0 - @eSquared * sin_pow_2(@phiPrime)) ** -1.5)
+       @etaSquared = (@v / @rho) - 1.0
+       @VII = Math.tan(@phiPrime) / (2 * @rho * @v)
+       @VIII =
+         (Math.tan(@phiPrime) / (24.0 * @rho * (@v ** 3.0)))\
+           * (5.0\
+             + (3.0 * tan_pow_2(@phiPrime))\
+             + @etaSquared\
+             - (9.0 * tan_pow_2(@phiPrime) * @etaSquared))
+       @IX =
+         (Math.tan(@phiPrime) / (720.0 * @rho * (@v ** 5.0)))\
+           * (61.0\
+             + (90.0 * tan_pow_2(@phiPrime))\
+             + (45.0 * tan_pow_2(@phiPrime) * tan_pow_2(@phiPrime)))
+       @X = sec(@phiPrime) / @v
+       @XI =
+         (sec(@phiPrime) / (6.0 * @v * @v * @v))\
+           * ((@v / @rho) + (2 * tan_pow_2(@phiPrime)))
+       @XII =
+         (sec(@phiPrime) / (120.0 * (@v ** 5.0)))\
+           * (5.0\
+             + (28.0 * tan_pow_2(@phiPrime))\
+             + (24.0 * tan_pow_2(@phiPrime) * tan_pow_2(@phiPrime)))
+       @XIIA =
+         (sec(@phiPrime) / (5040.0 * (@v ** 7.0)))\
+           * (61.0\
+             + (662.0 * tan_pow_2(@phiPrime))\
+             + (1320.0 * tan_pow_2(@phiPrime) * tan_pow_2(@phiPrime))\
+             + (720.0\
+               * tan_pow_2(@phiPrime)\
+               * tan_pow_2(@phiPrime)\
+               * tan_pow_2(@phiPrime)))
+       @phi =
+         @phiPrime\
+           - (@VII * ((@E - @E0) ** 2.0))\
+           + (@VIII * ((@E - @E0) ** 4.0))\
+           - (@IX * ((@E - @E0) ** 6.0))
+       @lambda =
+         @lambda0\
+           + (@X * (@E - @E0))\
+           - (@XI * ((@E - @E0) ** 3.0))\
+           + (@XII * ((@E - @E0) ** 5.0))\
+           - (@XIIA * ((@E - @E0) ** 7.0))
+
+      { :latitude => rad_to_deg(@phi), :longitude => rad_to_deg(@lambda) }
+       
+    end
+    
+    # Takes OSGB36 Latitude and Longitude coords 
+    # and returns WGS84 Latitude and Longitude
+    def to_WGS84(latitude,longitude)
+      
+      @a        = 6377563.396
+      @b        = 6356256.909
+      @eSquared = ((@a * @a) - (@b * @b)) / (@a * @a)
+
+      @phi = deg_to_rad(latitude)
+      @lambda = deg_to_rad(longitude)
+      @v = @a / (Math.sqrt(1 - @eSquared * sin_pow_2(@phi)))
+      @H = 0
+      @x = (@v + @H) * Math.cos(@phi) * Math.cos(@lambda)
+      @y = (@v + @H) * Math.cos(@phi) * Math.sin(@lambda)
+      @z = ((1 - @eSquared) * @v + @H) * Math.sin(@phi)
+
+      @tx =        446.448
+      @ty =       -124.157
+      @tz =        542.060
+      
+      @s  =         -0.0000204894
+      @rx = deg_to_rad( 0.00004172222)
+      @ry = deg_to_rad( 0.00006861111)
+      @rz = deg_to_rad( 0.00023391666)
+
+      @xB = @tx + (@x * (1 + @s)) + (-@rx * @y)     + (@ry * @z)
+      @yB = @ty + (@rz * @x)      + (@y * (1 + @s)) + (-@rx * @z)
+      @zB = @tz + (-@ry * @x)     + (@rx * @y)      + (@z * (1 + @s))
+
+      @a        = 6378137.000
+      @b        = 6356752.3141
+      @eSquared = ((@a * @a) - (@b * @b)) / (@a * @a)
+
+      @lambdaB = rad_to_deg(Math.atan(@yB / @xB))
+      @p = Math.sqrt((@xB * @xB) + (@yB * @yB))
+      @phiN = Math.atan(@zB / (@p * (1 - @eSquared)))
+      
+      (1..10).each do |i|
+        @v = @a / (Math.sqrt(1 - @eSquared * sin_pow_2(@phiN)))
+        @phiN1 = Math.atan((@zB + (@eSquared * @v * Math.sin(@phiN))) / @p)
+        @phiN = @phiN1
+      end
+
+      @phiB = rad_to_deg(@phiN)
+
+      { :latitude => @phiB, :longitude => @lambdaB }
+         
+    end
+    
+    private # Some Math
+    
+    def deg_to_rad(degrees)
+      degrees / 180.0 * Math::PI
+    end
+
+    def rad_to_deg(r)
+      (r/Math::PI)*180
+    end
+
+    def sin_pow_2(x)
+      Math.sin(x) * Math.sin(x)
+    end
+
+    def cos_pow_2(x)
+      Math.cos(x) * Math.cos(x)
+    end
+
+    def tan_pow_2(x)
+      Math.tan(x) * Math.tan(x)
+    end
+
+    def sec(x)
+      1.0 / Math.cos(x)
+    end
+  
+end

scripts/en2latlon.rb

@@ -0,0 +1,31 @@
+# Convert OSGB36 eastings and northings to WGS84 latitudes and longitudes
+# $ ruby en2latlon.rb infile.csv > outfile.csv
+# http://adrianshort.org/2013/01/23/converting-eastings-and-northings-to-latitudes-and-longitudes-in-ruby/
+require 'csv'
+require_relative './OSGB36'
+
+lines = CSV.read(ARGV[0]) # read the whole file into an array of arrays
+
+easting_col = nil
+northing_col = nil
+
+(0..lines.size - 1).each do |i|
+  if i == 0
+    easting_col = lines[i].index("Eastings")
+    northing_col = lines[i].index("Northings")
+    lines[i] << "lat"
+    lines[i] << "lng"
+  else
+    ll = OSGB36.en_to_ll(lines[i][easting_col].to_f, lines[i][northing_col].to_f)
+    lines[i] << ll[:latitude]
+    lines[i] << ll[:longitude]
+  end
+end
+
+csv_str = CSV.generate do |csv|
+  lines.each do |line|
+    csv << line
+  end
+end
+
+puts csv_str

scripts/postcode-import.rb

@@ -0,0 +1,32 @@
+# Import postcode data from CodePoint Open
+
+require 'csv'
+require_relative "../app"
+
+i = 0
+
+CSV.foreach(ARGV.shift) do |row|
+  i += 1
+  next if i == 1 # skip header row
+  puts i
+  puts row
+  puts 
+  
+  ward = District.first(:ons_district_code => row[9])
+  
+  puts Postcode.create!(
+    :postcode => row[0],
+    :positional_quality_indicator => row[1],
+    :eastings => row[2],
+    :northings => row[3],           
+    :country_code => row[4],          
+    :nhs_regional_ha_code => row[5],      
+    :nhs_ha_code => row[6],
+    :admin_county_code => row[7],        
+    :admin_district_code => row[8],  
+    :admin_ward_code => row[9],
+    :lat => row[10],    
+    :lng => row[11],               
+    :ward_id => ward.id
+  )    
+end

scripts/wardcodes.rb

@@ -0,0 +1,14 @@
+# Add new ONS ward codes to districts table
+
+# E05000555
+
+require_relative "../app"
+
+c = 5000555
+
+District.all(:body_id => 1, :order => [:name]).each do |d|
+  d.ons_district_code = "E%08d" % c
+  d.save
+  puts d.name, d.ons_district_code
+  c += 1
+end