UK Postcode distance 4

[spperl]

Hi,

I've been on websites where it allows you to filter results based on a specified distance from a provided postcode.

Two questions I have are as follows:

1. What is needed to make this possible, and,

2. Are there any perl modules out there that can assist?

Any help would be gratefully appreciated and in this instance may lead to potential work.

I look forward to your replies.

Thanks you.

 

[Kirsle]

I have working code for this already:

Demo:

http://www.cuvou.com/wizards/distance.cgi

Source:

http://www.cuvou.com/wizards/distance.pls

Zipcode database:

http://www.cuvou.com/wizards/ZIP_CODES.txt

Zipcodes as Perl structure:

http://www.cuvou.com/wizards/zip.txt

The interesting parts of the CGI:

# Calculate the distance between two points.
sub distance {
	my ($lat1, $lon1, $lat2, $lon2, $unit) = @_;
	my $theta = $lon1 - $lon2;
	my $dist = sin(deg2rad($lat1)) * sin(deg2rad($lat2)) + cos(deg2rad($lat1)) * cos(deg2rad($lat2)) * cos(deg2rad($theta));
	$dist  = acos($dist);
	$dist = rad2deg($dist);
	$dist = $dist * 60 * 1.1515;
	if ($unit eq "K") {
		$dist = $dist * 1.609344;
	}
	elsif ($unit eq "N") {
		$dist = $dist * 0.8684;
	}
	return ($dist);
}

# This function get the arccos function using arctan function
sub acos {
	my ($rad) = @_;
	my $ret = atan2(sqrt(1 - $rad**2), $rad);
	return $ret;
}


# This function converts decimal degrees to radians
sub deg2rad {
	my ($deg) = @_;
	return ($deg * $pi / 180);
}

# This function converts radians to decimal degrees
sub rad2deg {
	my ($rad) = @_;
	return ($rad * 180 / $pi);
}

# Round numbers down.
sub round {
	my $num = shift;

	if ($num =~ /\./) {
		$num =~ s~^(\d+)\.(\d\d)\d+$~$1.$2~ig;
	}

	return $num;
}

The first function there gives ya the distance between two latitude/longitude coordinates. The zipcode DB just provides the coordinates for every zipcode in the U.S.

I got the DB file from a site that sells a commercial product for doing exactly this... they gave out their zipcode database for free along with those Perl functions. The DB might be a little bit outdated; if you find a newer DB that gives coordinates for all zip codes, feel free to use that instead.

Cuvou.com | My personal homepage

perl -e '$|=$i=1;print" oo\n<|>\n_|_";x:sleep$|;print"\b",$i++%2?"/":"_";goto x;'

 

[Kirsle]

*Ah, I read "UK" as "US"

The code is the same, you'll just need to find a database of UK postcodes and their coordinates.

Cuvou.com | My personal homepage

perl -e '$|=$i=1;print" oo\n<|>\n_|_";x:sleep$|;print"\b",$i++%2?"/":"_";goto x;'

 

[spperl]

This was much, much more than I was anticipating so thanks so much.

 

[ChrisHunt]

What is needed to make this possible

You'll need some way of getting from a postcode to geographical coordinates. That's not as easy as it should be, as the Post Office charge a lot of cash to provide that service, even though we've already paid for it with our taxes (more on this issue at

http://www.freeourdata.org.uk

).

You can get postcode coordinates for free from

http://ernestmarples.com

, they also have a perl routine at

http://search.cpan.org/dist/WWW-ErnestMarples/

Once you've got a lat,long pair, you have to determine the distance between. Kirsle's procedure looks suspiciously simple. I'd want to be sure that it applies everywhere in the world instead of just the US.

The approach I'd use would be to use the Geography::NationalGrid library to convert from lat/long to an OS map reference. That will give you an x,y coordinate expressed in kilometers, from which you can use simple trigonometry to work out the distance.

However, if you want 100% accuracy you may need to convert between WGS-84 lat/long (which is the worldwide standard used for GPS, online maps etc) and the Airey/OSGB36 standard used by the Ordnance survey. You can read about that at

http://www.dorcus.co.uk/carabus/ll_ngr.html

. I based the following perl code on the javascript on that page:

sub ll_to_cartesian {
   my $lat = shift;
   my $lon = shift;
   my $axis = shift;
   my $ecc = shift;
   my $height = shift;

   my $v = $axis / (sqrt(1 - $ecc * (sin($lat)**2)));
   my $x = ($v + $height) * cos($lat) * cos($lon);
   my $y = ($v + $height) * cos($lat) * sin($lon);
   my $z = ((1 - $ecc) * $v + $height) * sin($lat);
   return ($x,$y,$z);
}

sub cartesian_to_ll {
   my $x = shift;
   my $y = shift;
   my $z = shift;
   my $axis = shift;
   my $ecc = shift;

   my $lon = atan2($y,$x);
   my $p = sqrt(($x * $x) + ($y * $y));
   my $lat = atan2($z,($p * (1 - $ecc)));
   my $v = $axis / (sqrt(1 - $ecc * (sin($lat) * sin($lat))));
   my $errvalue = 1.0;
   my $lat0 = 0;
   while ($errvalue > 0.001) {
      $lat0 = atan2(($z + $ecc * $v * sin($lat)), $p);  
      $errvalue = abs($lat0 - $lat);
      $lat=$lat0;
   }
   my $height = $p / cos($lat) - $v;
   return($lat,$lon);
}

sub airey_to_wgs84 {
   my $latin = shift;
   my $longin = shift;

   $latin = Geography::NationalGrid->deg2rad($latin);
   $longin = Geography::NationalGrid->deg2rad($longin);

   my ($x,$y,$z) = ll_to_cartesian ($latin,$longin,6377563.396,0.00667054,0);

   $x += 371;
   $y -= 112;
   $z += 434;

   my ($lat,$long) = cartesian_to_ll ($x,$y,$z,6378137.0,0.00669438);

   $lat = Geography::NationalGrid->rad2deg($lat);
   $long = Geography::NationalGrid->rad2deg($long);

   return ($lat,$long);
}

sub wgs84_to_airey_ {
   my $latin = shift;
   my $longin = shift;

   $latin = Geography::NationalGrid->deg2rad($latin);
   $longin = Geography::NationalGrid->deg2rad($longin);

   my ($x,$y,$z) = ll_to_cartesian ($latin,$longin,6378137.0,0.00669438,0);

   $x -= 371;
   $y += 112;
   $z -= 434;

   my ($lat,$long) = cartesian_to_ll ($x,$y,$z,6377563.396,0.00667054);

   $lat = Geography::NationalGrid->rad2deg($lat);
   $long = Geography::NationalGrid->rad2deg($long);

   return ($lat,$long);
}

So you get a lat/long from Ernest Maples, feed it into wgs_to_airey(), and feed the result into Geography::NationalGrid to get an easting/northing reference.

Maybe you're better with Kirsle's code...

-- Chris Hunt
Webmaster & Tragedian
Extra Connections Ltd