Skip to content

Instantly share code, notes, and snippets.

@turugina
Last active August 29, 2015 14:22
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save turugina/eda38cd9a10a5c0bc28f to your computer and use it in GitHub Desktop.
Save turugina/eda38cd9a10a5c0bc28f to your computer and use it in GitHub Desktop.
緯度経度3点から大体の面積を求める
use strict;
use warnings;
use utf8;
use Math::Trig;
use v5.12;
# lat 緯度 lng 経度
my $p1 = [ 34.678506,135.498941 ];
my $p2 = [ 34.6785, 135.499194 ];
my $p3 = [ 34.677377,135.499024 ];
# 長半径[m]
my $A = 6378137;
# 扁平率
#my $F = 1.0 / 298.257222101;
# 円周[m]
my $L = 2 * 3.14159 * $A;
say calcs($p1, $p2, $p3);
say calcs($p2, $p3, $p1);
say calcs($p3, $p1, $p2);
sub calcunit {
my ($lat, $lng) = @_;
my $ulat = $L/360.0;
# 北緯xx度地点での半径
my $r = $A * cos(deg2rad($lat));
#say "r($lat) = $r";
my $l = 2 * 3.14159 * $r;
my $ulng = $l/360.0;
return ($ulat, $ulng);
}
sub calcs {
my ($p1, $p2, $p3) = @_;
my ($ulat,$ulng) = calcunit(($p1->[0]+$p2->[0]+$p3->[0])/3.0, ($p1->[1]+$p2->[1]+$p3->[1])/3.0);
#printf("p1 unit lat m/deg = %.2f\n", $ulat);
#printf("p1 unit lng m/deg = %.2f\n", $ulng);
#printf("p1-p2 lat diff = %f\n", abs($p2->[0]-$p1->[0]));
#printf("p1-p2 lng diff = %f\n", abs($p2->[1]-$p1->[1]));
my $P2 = [($p2->[0]-$p1->[0])*$ulat, ($p2->[1]-$p1->[1])*$ulng];
my $P3 = [($p3->[0]-$p1->[0])*$ulat, ($p3->[1]-$p1->[1])*$ulng];
#say join(",", map{int} @$P2);
#say join(",", map{int} @$P3);
0.5 * abs($P2->[0]*$P3->[1] - $P3->[0]*$P2->[1]);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment