OsmocomBB

#!/usr/bin/perl

#############################################################################

## >>> 2014 <<<<

## Martin Hauke (mardnh@gmx.de)

##

## Older GSM BTS's from O2 Germany (formerly MNC 07 - now MNC 03) send out their BTS

## coordinates via smscb (channel 221).

## With Alex Badea's patches for libosmocore/osmocombb applied you can log those smscb

## messages via osmocom-bb's cell_log application.

## This script parses osmocom-bb cell_log output and writes the relevant info to a

## nice kml-file.

##

## usage example:

## cat cell.log | o2-osmocom-bb-smscb2kml.pl > nice_map.kml

##

##

## The used Converter for Gauss - Krueger coordinates is licensed under the GPL and was

## written by

## - Andreas Achtzehn (achtzehn@linux-society.de)

## - Nobert Hüttisch (nobbi@nobbi.com)

#use strict;

my %seen;

my ($line, $time, $arfcn, $mcc, $mnc, $lac, $cid, $cb);

my ($viag_gk, $sr, $sx);

my $regex = qw(^CellBroadcast: TIME=(?<time>\d+)\sARFCN=(?<arfcn>\d+)\sMCC=(?<mcc>\d+)\sMNC=(?<mnc>\d+)\sLAC=(?<lac>.*)\sCID=(?<cid>.*)\sCB=\"(?<cb>.*)\"$);

&printKMLheader;

foreach $line (<STDIN>) {

	chomp($line);

	if ( $line =~ /$regex/ ) {

		#print "$line\n";

		$time  = "$+{time}";

		$arfcn = "$+{arfcn}";

		$mcc   = "$+{mcc}";

		$mnc   = "$+{mnc}";

		$lac   = "$+{lac}";

		$cid   = "$+{cid}";

		$cb    = "$+{cb}";

		$cid_dec = hex($cid);

		$lac_dec = hex($lac);

		$time_x  = scalar localtime($time);

		# check if O2-Germany (O2 MNC was until mid 2016 - since then 03)

		if ($mnc == "03" && $mcc == "262") {

			# did we seen that entry point already?

			next if $seen{"$arfcn,$mcc,$mnc,$lac,$cid"}++;

			$cb =~ /(?<gk>\d+).*/;

			$viag_gk = "$+{gk}";

			$xlat=substr($viag_gk,0,6)*10;

			$xlon=substr($viag_gk,6,6)*10;

			getCoordinate($xlat, $xlon);

			#print "Google Maps: https://maps.google.com/maps?z=12&h=m&q=loc:$wgs84lat+$wgs84lon&z=15\n";

			#print "$time $arfcn $mcc $mnc $lac $cid $viag_gk $wgs84lat $wgs84lon\n";

			# print placemarks

			print qq{\t\t<Placemark>\n};

			print qq{\t\t\t<name>CID:$cid_dec</name>\n};

			print qq{\t\t\t<description><![CDATA[CID: <b>$cid ($cid_dec)</b><br> LAC: <b>$lac ($lac_dec)</b><br>ARFCN: <b>$arfcn</b><br> MCC: <b>$mcc</b><br>MNC: <b>$mnc</b><br>TIME: <b>$time_x</b> ]]></description>\n};

			#print qq{\t\t\t</description>\n};

			#print "<styleUrl>#rxl69</styleUrl>"

			print qq{\t\t\t<Point>\n};

			print qq{\t\t\t\t<coordinates>$wgs84lon,$wgs84lat</coordinates>\n};

			print qq{\t\t\t</Point>\n};

			print qq{\t\t</Placemark>\n};

		};

	};

};

&printKMLfooter;

sub getCoordinate {

	$sr = shift;

	$sx = shift;

	&potsdam;

	&modolensky;

}

sub potsdam {

	# Errechnen der Position nach Potsdam-Koordinatensystem

	$bm  = int($sr/1000000);

	$y   = $sr-($bm*1000000+500000);

	$si  = $sx/111120.6196;

	$px  = $si+0.143885358*sin(2*$si*0.017453292)+0.00021079*sin(4*$si*0.017453292)+0.000000423*sin(6*$si*0.017453292);

	$t   = (sin($px*0.017453292))/(cos($px*0.017453292));

	$v   = sqrt(1+0.006719219*cos($px*0.017453292)*cos($px*0.017453292));

	$ys  = ($y*$v)/6398786.85;

	$lat = $px-$ys*$ys*57.29577*$t*$v*$v*(0.5-$ys*$ys*(4.97-3*$t*$t)/24);

	$dl  = $ys*57.29577/cos($px*0.017453292) * (1-$ys*$ys/6*($v*$v+2*$t*$t-$ys*$ys*(0.6+1.1*$t*$t)*(0.6+1.1*$t*$t)));

	$lon = $bm*3+$dl;

	$ph =int($lat);

	$pz =($lat-$ph)*60;

	$pm =int($pz);

	$ps =($pz-$pm)*60;

	$lh =int($lon);

	$lz =($lon-$lh)*60;

	$lm =int($lz);

	$ls =($lz-$lm)*60;

}

sub modolensky {

	$potsd_a  = 6377397.155;

	$wgs84_a  = 6378137.0;

	$potsd_f  = 1/299.152812838;

	$wgs84_f  = 1/298.257223563;

	$potsd_es = 2*$potsd_f - $potsd_f*$potsd_f;

	$potsd_dx = 606.0;

	$potsd_dy = 23.0;

	$potsd_dz = 413.0;

	$pi = 3.141592654;

	$latr = $lat/180*$pi;

	$lonr = $lon/180*$pi;

	$sa = sin($latr);

	$ca = cos($latr);

	$so = sin($lonr);

	$co = cos($lonr);

	$bda  = 1-$potsd_f;

	$delta_a = $wgs84_a - $potsd_a;

	$delta_f = $wgs84_f - $potsd_f;

	$rn = $potsd_a / sqrt(1-$potsd_es*sin($latr)*sin($latr));

	$rm = $potsd_a * ((1-$potsd_es)/sqrt(1-$potsd_es*sin($latr)*sin($latr)*1-$potsd_es*sin($latr)*sin($latr)*1-$potsd_es*sin($latr)*sin($latr)));

	$ta = (-$potsd_dx*$sa*$co - $potsd_dy*$sa*$so)+$potsd_dz*$ca;

	$tb = $delta_a*(($rn*$potsd_es*$sa*$ca)/$potsd_a);

	$tc = $delta_f*($rm/$bda+$rn*$bda)*$sa*$ca;

	$dlat = ($ta+$tb+$tc)/$rm;

	$dlon = (-$potsd_dx*$so + $potsd_dy*$co)/($rn*$ca);

	$wgs84lat = ($latr + $dlat)*180/$pi;

	$wgs84lon = ($lonr + $dlon)*180/$pi;

	$ph=int($wgs84lat);

	$pz=($wgs84lat-$ph)*60;

	$pm=int($pz);

	$ps=int(($pz-$pm)*60);

	$lh=int($wgs84lon);

	$lz=($wgs84lon-$lh)*60;

	$lm=int($lz);

	$ls=int(($lz-$lm)*60);

	return ($wgs84lat, $wgs84lon )

}

sub printKMLheader {

	print qq{<?xml version="1.0" encoding="UTF-8"?>\n};

	print qq{<kml xmlns="http://www.opengis.net/kml/2.2">\n};

	print qq{\t<Document>\n};

};

sub printKMLfooter {

	print qq{\t</Document>\n};

	print qq{</kml>\n};

};

# EOF