(Library) 緯度・経度情報から面積を計算する

作成者: 西本絵梨子

編集履歴:

• 2018-8-24(Fri): 新規作成

緯度、経度の情報から面積(≒矩形面積)を計算します。

緯度幅を小さく取ることで、矩形近似をしています。

間違いなどを見つけられましたら、西本( eriko _at_ gfd-dennou.org )までご連絡いただけると助かります。

calc_area_from_lat_lon.rb

• GPhys

経度slonからelon, 緯度slatからelatで囲まれた部分の面積を計算。

GPhys_Function::calc_area_from_lat_lon( slon, elon, slat, elat, n_dlat )

ARGUMENTS
 * slon, elon   : start and end points in longitude (Numeric or UNumeric)
 * slat, elat   : start and end points in latitude (Numeric or UNumeric)
 * n_dlat       : number of partitions in latitude (Numeric, default=10)

RETURN VALUE
 * area         : area enclosed by slon, elon, slat, elat (UNumeric)

require "numru/gphys"
include NumRu

module GPhys_Function

  module_function

  #-- constant
  @radius=UNumeric.new(6.37E6,"m") # radius of the planet

  def deg_to_rad( deg, type )
    if deg.class!=UNumeric
      if type=="lon"
        deg=UNumeric.new(deg,"degrees_east")
      elsif type=="lat"
        deg=UNumeric.new(deg,"degrees_north")
      end
    end
    rad=deg.convert("rad")

    return rad, deg
  end#-- def deg_to_rad( deg, type )

  def calc_area_from_lat_lon( slon, elon, slat, elat, n_dlat=10 )
    slon_rad, = deg_to_rad( slon, "lon" )
    elon_rad, = deg_to_rad( elon, "lon" )
    slat_rad, = deg_to_rad( slat, "lat" )
    elat_rad, = deg_to_rad( elat, "lat" )

    #-- dlon
    lon_ratio=(elon_rad-slon_rad)/(UNumeric.new(2.0*Math::PI,"rad"))

    #-- Δlatを小さくとって緯度円計算し、足し合わせ
    dlat_rad=(elat_rad-slat_rad)/n_dlat.to_f
    lat_length=@radius*dlat_rad.val

    puts "n_dlat=#{n_dlat}\n"
    puts "dlat=#{dlat_rad}, #{dlat_rad.convert("degrees_north")}\n"
    puts "lat_length=#{lat_length}\n\n"

    area=0.0
    n_dlat.times{|ilat|
      theta=slat_rad+dlat_rad*ilat
      lon_circ_length=(2.0*Math::PI*@radius*theta.cos)*lon_ratio

      area_tmp=lon_circ_length*lat_length
#p "#{theta.convert("degrees_north")}, #{lon_circ_length}, area=#{area_tmp}"
      if ilat==0
        area=area_tmp
      else
        area+=area_tmp
      end
    }#-- n_dlat.times{|ilat|

    return area

  end#-- def calc_area_in_lat_lon( slon, elon, slat, elat, n_dlat=nil )


  #----- Sample
  if __FILE__==$0

    slon, elon, slat, elat = 0, 360, -90, 90
    n_dlat=180
    area=calc_area_from_lat_lon( slon, elon, slat, elat, n_dlat )

    areaEarth=4.0*Math::PI*@radius**2 

    puts "\nerror percentage"
    puts "#{((area-areaEarth)/areaEarth)*100} %"

  end

end#-- module GPhys_Function