Convert RT90 to WGS84

geo_conv.rb

def convert_rt90_to_wgs84(geopos)
  # from http://mellifica.se/geodesi/gausskruger.js
  geopos=geopos.first

  x = geopos['x'].to_i
  y = geopos['y'].to_i

  axis = 6378137.0 # GRS 80.
  flattening = 1.0 / 298.257222101 # GRS 80.
  central_meridian = 15.0 + 48.0/60.0 + 22.624306/3600.0
  scale = 1.00000561024
  false_northing = -667.711
  false_easting = 1500064.274


  lat_lon = [nil,nil]
  if (central_meridian.nil?)
    return lat_lon
  end
  # Prepare ellipsoid-based stuff.
  e2 = flattening * (2.0 - flattening)
  n = flattening / (2.0 - flattening)
  a_roof = axis / (1.0 + n) * (1.0 + n*n/4.0 + n*n*n*n/64.0)
  delta1 = n/2.0 - 2.0*n*n/3.0 + 37.0*n*n*n/96.0 - n*n*n*n/360.0
  delta2 = n*n/48.0 + n*n*n/15.0 - 437.0*n*n*n*n/1440.0
  delta3 = 17.0*n*n*n/480.0 - 37*n*n*n*n/840.0
  delta4 = 4397.0*n*n*n*n/161280.0

  astar = e2 + e2*e2 + e2*e2*e2 + e2*e2*e2*e2
  bstar = -(7.0*e2*e2 + 17.0*e2*e2*e2 + 30.0*e2*e2*e2*e2) / 6.0
  cstar = (224.0*e2*e2*e2 + 889.0*e2*e2*e2*e2) / 120.0
  dstar = -(4279.0*e2*e2*e2*e2) / 1260.0

  # Convert.
  deg_to_rad = Math::PI / 180
  lambda_zero = central_meridian * deg_to_rad
  xi = (x - false_northing) / (scale * a_roof)
  eta = (y - false_easting) / (scale * a_roof)
  xi_prim = xi -
      delta1*Math.sin(2.0*xi) * Math.cosh(2.0*eta) -
      delta2*Math.sin(4.0*xi) * Math.cosh(4.0*eta) -
      delta3*Math.sin(6.0*xi) * Math.cosh(6.0*eta) -
      delta4*Math.sin(8.0*xi) * Math.cosh(8.0*eta)
  eta_prim = eta -
      delta1*Math.cos(2.0*xi) * Math.sinh(2.0*eta) -
      delta2*Math.cos(4.0*xi) * Math.sinh(4.0*eta) -
      delta3*Math.cos(6.0*xi) * Math.sinh(6.0*eta) -
      delta4*Math.cos(8.0*xi) * Math.sinh(8.0*eta)
  phi_star = Math.asin(Math.sin(xi_prim) / Math.cosh(eta_prim))
  delta_lambda = Math.atan(Math.sinh(eta_prim) / Math.cos(xi_prim))
  lon_radian = lambda_zero + delta_lambda
  lat_radian = phi_star + Math.sin(phi_star) * Math.cos(phi_star) *
      (astar +
          bstar*(Math.sin(phi_star) ** 2) +
          cstar*(Math.sin(phi_star) ** 4) +
          dstar*(Math.sin(phi_star) ** 6))
  lat_lon[0] = ((lat_radian * 1800000.0 / Math::PI).round).to_f/10000
  lat_lon[1] = ((lon_radian * 1800000.0 / Math::PI).round).to_f/10000
  return lat_lon
end