File Coverage

blib/lib/Geo/Coordinates/Convert.pm

Criterion	Covered	Total	%
statement	38	38	100.0
branch	1	2	50.0
condition			n/a
subroutine	8	8	100.0
pod	0	3	0.0
total	47	51	92.1

package Geo::Coordinates::Convert;

use 5.006;
use strict;
use warnings;

use Exporter ();
use vars qw($VERSION @ISA @EXPORT);

@EXPORT = qw(
	set_Mean_Longitude
	geo2lII
	lII2geo
);

@ISA = qw(Exporter);

$VERSION = '0.01';

# some constants
my $PIsur4 = atan2(1, 1);
my $PIsur180 = $PIsur4 / 45;
my $PI = $PIsur4 * 4;
my $PIsur2 = $PIsur4 * 2;
my $precision = 0.001; # degrees

# these "constants" are relative to France (see www.ign.fr)
my $N = 0.7289686274;
my $C = 11745793.39;
my $XS = 600000;
my $YS = 8199695.768; # Lambert II (6199695.768 for extended LambertII)
my $E = 0.08248325676;
my $Esur2 = $E / 2;

# mean longitude. 2.3372 is the longitude of the meridian of Paris (degrees)
my $mean_longitude = 2.3372;

sub set_Mean_Longitude {
	$mean_longitude = $_[0];
}

sub geo2lII {
# Conversion from geographics coordinates (decimal degrees) to Lambert II coordinates (meters)
#
# list in (decimal degrees):
#   (geographic longitude, geographic latitude, mean longitude)
# list out (meters):
#   (lambert II longitude, lambert II latitude)
#
# longitudes from -180° west to +180° east
# latitudes from -90° (south) to +90° (north)

my ($l, $r, $gamma, $sinlat, $xl, $yl);
	# arguments
	my ($lon, $lat, $ml) = @_;
	$sinlat = sin($lat * $PIsur180);
	set_Mean_Latitude($ml) if $ml;

	$l = 0.5 * (
		log (
			(1 + $sinlat) / (1 - $sinlat)
		)
		- $E * log (
			(1 + $E * $sinlat) / (1 - $E * $sinlat)
		)
	);

	$r = $C * exp(-$N * $l);

	$gamma = $N * ($lon - $mean_longitude) * $PIsur180;

	$xl = $XS + $r * sin($gamma);
	$yl = $YS - $r * cos($gamma);

	return ($xl, $yl);
}

sub lII2geo {
# Conversion from Lambert II coordinates (meters) to geographics coordinates (decimal degrees)
#
# list in (meters):
#   (lambert II longitude, lambert II latitude)
# list out (decimal degrees):
#   (geographic longitude, geographic latitude)
# conventions:
# - negatives western longitudes, positives eastern
# - negatives southern latitudes, positives northern

	my ($el, $r, $gamma, $lon, $lat, $latp);

	my ($xl, $yl) = @_;

	$r = sqrt( ($xl - $XS) * ($xl - $XS) + ($yl - $YS) * ($yl - $YS) );
	$gamma = atan2 ( ($xl - $XS) / ($YS - $yl), 1 );

	$lon = $mean_longitude * $PIsur180 + $gamma / $N;
	$el = exp( - log ($r / $C) / $N );

	$lat = $PIsur4;
	$latp = 0;
	while (abs($latp - $lat) > $precision) {
		$latp = $lat;
		$lat = 2 * atan2( ((1 + $E * $latp) / (1 - $E * $latp) )** $Esur2 * $el, 1 ) - $PIsur2;
	}

	return ($lon/$PIsur180, $lat/$PIsur180);
}


1;
__END__