Inverse Distance Weighted (IDW) Interpolation on Earth Surface in Ruby

earth_interpolation.rb

class Weather::Interpolation
  def initialize(x, y, z, xi, yi)
    @x, @y, @z, @xi, @yi = x, y, z, xi, yi
  end

  def calc
    # if any location closer than 10km, return this location
    close_location = distances.index {|d| d < 10}
    close_location ? @z[close_location] : interpolation
  end

  private
  def interpolation
    weights.each_with_index.map do |wi, i|
      val = @z[i]
      val = @z.compact.sum / @z.compact.size if val.nil? # use mean if nil
      wi * val
    end.sum
  end

  def weights
    v = distances.map { |d| 1 / (d**2) }
    sum = v.sum
    v.map { |vi| vi / sum }
  end

  def distances
    @distances ||= @x.each_with_index.map do |lat, i|
      lng = @y[i]
      distance(lat, lng, @xi, @yi)
    end
  end

  # https://stackoverflow.com/questions/12966638/
  # how-to-calculate-the-distance-between-two-gps-coordinates-without-using-google-m
  def distance(lat1, lng1, lat2, lng2)
    rad_per_deg = Math::PI/180  # PI / 180
    rkm = 6371                  # Earth radius in kilometers
    rm = rkm * 1000             # Radius in meters

    dlat_rad = (lat2-lat1) * rad_per_deg  # Delta, converted to rad
    dlon_rad = (lng2-lng1) * rad_per_deg

    lat1_rad, lon1_rad = lat1 * rad_per_deg, lng1 * rad_per_deg
    lat2_rad, lon2_rad = lat2 * rad_per_deg, lng2 * rad_per_deg

    a = Math.sin(dlat_rad/2)**2 + Math.cos(lat1_rad) * Math.cos(lat2_rad) * Math.sin(dlon_rad/2)**2
    c = 2 * Math::atan2(Math::sqrt(a), Math::sqrt(1-a))

    rm * c / 1000 # Delta in kilometers
  end
end