Import original source of Map-Splat 0.01

1 files changed, 521 insertions, 0 deletions

diff --git a/Splat.pm b/Splat.pm
new file mode 100644
index 0000000..5729345
--- /dev/null
+++ b/Splat.pm
@@ -0,0 +1,521 @@
+package Map::Splat;
+
+use 5.006;
+use strict;
+use warnings FATAL => 'all';
+
+use File::chdir;
+use File::Temp qw(tempdir);
+use File::Path qw(make_path);
+use File::Slurp qw(read_file write_file);
+use XML::LibXML::Simple qw(XMLin);
+use List::Util qw(min max);
+use Data::Dumper;
+use Image::Magick;
+use Moo;
+use POSIX ('floor', 'ceil');
+use v5.10;
+
+=head1 NAME
+
+Map::Splat - Wrapper for the Splat! RF propagation modeling tool
+
+=head1 VERSION
+
+Version 0.01
+
+=cut
+
+our $VERSION = '0.01';
+
+=head1 SYNOPSIS
+
+    use Map::Splat;
+
+    my $foo = Map::Splat->new();
+    ...
+
+=head1 METHODS
+
+=over 4
+
+=cut
+
+# requires:
+# - Splat
+# - wget
+# - zip
+# - gdal-bin package (not used yet, but just wait)
+# - XML::LibXML
+# - ImageMagick
+
+my $SRTM_DIR = $ENV{HOME} . '/.splat';
+my $SRTM_PATH = 'http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America';
+my $SPLAT_CMD = '/usr/bin/splat';
+
+=item new PARAMETERS
+
+Constructs a new mapping scenario. This takes the following parameters:
+
+lon: the tower longitude (east is positive)
+
+lat: the tower latitude (north is positive)
+
+height: the tower height above ground level, in feet
+
+freq: the frequency, in MHz
+
+h_width: the horizontal beamwidth, in degrees
+
+v_width: the vertical beamwidth, in degrees
+
+azimuth: the antenna orientation, in degrees from north (90 = east)
+
+tilt: the downward tilt of the antenna, in degrees
+
+and optionally accepts:
+
+dir: the working directory to use; defaults to a random tempdir
+
+id: the filename prefix for working files; defaults to '1'
+
+name: a string describing the tower; defaults to 'Tower'
+
+min_loss: the lowest path loss level to show on the map; defaults to 80 dB
+
+max_loss: the highest path loss level to show on the map; defaults to 160 dB
+
+=cut
+
+foreach ( qw( lon lat height freq h_width v_width azimuth tilt ) ) {
+  has $_ => ( is => 'rw', required => 1 );
+}
+
+has 'dir' => (
+  is => 'ro',
+  default => sub {
+    tempdir(CLEANUP => 1)
+  },
+);
+
+has 'id' => (
+  is => 'ro',
+  default => sub { 1 }
+);
+
+has 'name' => (
+  is => 'ro',
+  default => sub { 'Tower' }
+);
+
+# there are more sophisticated ways to calculate this. for now we just expose
+# the ends of the scale as optional parameters.
+has 'min_loss' => (
+  is => 'ro',
+  default => sub { 80 }
+);
+
+has 'max_loss' => (
+  is => 'ro',
+  default => sub { 160 }
+);
+
+=item image
+
+After calling I<calculate> this will contain an L<Image::Magick> object with
+the raw image data.
+
+=cut
+
+has 'image' => (
+  is => 'rw',
+);
+
+=item box
+
+After calling I<calculate> this will contain a hashref of 'north', 'south',
+'east', and 'west' indicating the bounds of the map area.
+
+=cut
+
+has 'box' => (
+  is => 'rw',
+  default => sub { +{} }
+);
+
+=item db_colors
+
+Returns a hashref of (path loss level in dB) => (color level).
+
+=cut
+
+sub db_colors {
+  my $self = shift;
+  my %color;
+  # we have 32 color bands; find the width of a single band
+  my $db_step = ($self->max_loss - $self->min_loss) / 32;
+  for (my $i = 0; $i < 32; $i++) {
+    # + 1 here so that the highest band is _above_ the -db switch setting
+    my $db = ceil($self->max_loss - ($i * $db_step)) + 1;
+    my $color = 248 - ($i * 8);
+    $color{$db} = $color;
+  }
+  \%color;
+}
+
+sub ew {
+  my $lon = shift;
+  $lon < 0 ? 'W' : 'E';
+}
+
+sub ns {
+  my $lat = shift;
+  $lat < 0 ? 'S' : 'N';
+}
+
+sub debug {
+  my $self = shift;
+  my $arg = shift;
+  if (ref $arg) {
+    warn Dumper($arg);
+  } else {
+    warn "$arg\n";
+  }
+}
+
+=item calculate
+
+Invoke Splat to create a map. Currently, the map will cover the 1-degree
+quadrangle containing the tower location.
+
+=cut
+
+sub calculate {
+  my $self = shift;
+
+  my $signed_lat = $self->lat;
+  my $signed_lon = $self->lon;
+
+  my $id = $self->id;
+
+  local $CWD = $self->dir;
+
+  # write the site file
+  # note that Splat wants longitude in degrees west, from 0 to <360
+  # (this will also matter for SDF files)
+  my $splat_lon = -1 * $signed_lon;
+  $splat_lon += 360 if $splat_lon < 0;
+
+  $self->debug("creating QTH file");
+  write_file("$id.qth",
+    $self->name . "\n" .
+    $signed_lat . "\n" .
+    $splat_lon . "\n" .
+    $self->height . "\n"
+  );
+
+  # Longley-Rice parameters
+  my $freq = $self->freq;
+  my $pol = 1; # polarization = vertical
+  $self->debug("creating LRP file for $freq MHz");
+
+  my $lrp = <<EOF;
+15.000  ; Earth Dielectric Constant (Relative permittivity)
+0.005   ; Earth Conductivity (Siemens per meter)
+301.000 ; Atmospheric Bending Constant (N-Units)
+$freq   ; Frequency in MHz (20 MHz to 20 GHz)
+5       ; Radio Climate
+$pol    ; Polarization (0 = Horizontal, 1 = Vertical)
+0.50    ; Fraction of situations
+0.50    ; Fraction of time
+0   ; ERP
+EOF
+  write_file("$id.lrp", $lrp);
+
+  # Loss contour file: specifies the color mapping for the output
+  $self->debug("creating $id.lcf\n");
+  my $lcf = '';
+  my $colors = $self->db_colors;
+  for my $db (sort { $a <=> $b } keys %$colors) {
+    my $red = $colors->{$db};
+    $lcf .= "$db: $red, 0, 0\n";
+  }
+  write_file("$id.lcf", $lcf);
+
+  # AZ file: radiation pattern in horizontal plane
+  # should be able to take this as an input (and kinda can, by overriding
+  # antenna_pattern)
+  my $bearing = $self->azimuth;
+  my $downtilt = $self->tilt;
+  my ($az, $el) = $self->antenna_pattern;
+  $self->debug("creating $id.az for bearing $bearing\n");
+
+  open my $fh, '>', "$id.az";
+  printf $fh "%0.1f\n", ($bearing);
+  foreach my $t (sort { $a <=> $b } keys (%$az)) {
+    printf $fh "%d\t%0.6f\n", ($t, $az->{$t});
+  } 
+  
+  close $fh;
+      
+  # EL file: radiation pattern in vertical plane
+  $self->debug("creating $id.el for downtilt $downtilt\n");
+  open $fh, '>', "$id.el";
+  printf $fh "%0.1f\t%0.1f\n", ($downtilt, $bearing);
+  foreach my $t (sort { $a <=> $b } keys (%$el)) {
+    printf $fh "%0.1f\t%0.3f\n", ($t, $el->{$t}); 
+  } 
+  close $fh;
+
+  # make sure we have SRTM for the area around the antenna
+  # ("area around the antenna" should be configurable)
+  # (and support SRTM1/high def version? only works in North America)
+
+  for (my $dx = -2; $dx <= 2; $dx++) {
+    for (my $dy = -2; $dy <= 2; $dy++) {
+      my $splat_lon_east = floor($splat_lon) + $dx;
+      my $splat_lon_west = ($splat_lon_east + 1) % 360;
+      my $lat_south = floor($signed_lat) + $dy;
+      my $lat_north = $lat_south + 1;
+
+      next if $lat_north > 90 or $lat_south < -90; # haha
+
+      # the naming convention used by Splat's SDF files:
+      # south:north:east:west
+      # longitude is degrees west, from 0 to 359.
+      # note that the last one (from 1E to 0W) is 359:0, not 359:360
+      my $sdffile = "$SRTM_DIR/$lat_south:$lat_north:$splat_lon_east:$splat_lon_west.sdf";
+      if ( ! -f $sdffile ) {
+        if ( ! -d $SRTM_DIR ) {
+          make_path($SRTM_DIR); # dies on error
+        }
+        local $CWD = $SRTM_DIR;
+        # the USGS, though, just marks them as N/S and E/W. reference is the
+        # southwest corner.
+        # that's the negative of splat longitude, shifted to the range (-180, 180)
+        # or in the east, it's 360 - splat longitude
+        my $usgs_lon_west = -1 * $splat_lon_west;
+        $usgs_lon_west += 360 if $usgs_lon_west < -180;
+        my $hgtfile = ns($lat_south) . abs($lat_south) .
+                      ew($usgs_lon_west) . abs($usgs_lon_west) .
+                      '.hgt';
+        my $zipfile = $hgtfile . '.zip';
+        $self->debug("fetching $SRTM_PATH/$zipfile\n");
+        my $status = system('wget', "$SRTM_PATH/$zipfile");
+        if ( $status > 0 ) {
+          # normal if the tower location is near a coast, so try to continue
+          warn "couldn't download $zipfile\n";
+          next;
+        }
+        $status = system('unzip', $zipfile);
+        # this isn't normal though
+        die "couldn't extract $zipfile\n" if $status > 0;
+        # unfortunately gives no useful output
+        $status = system('srtm2sdf', $hgtfile);
+        unlink $hgtfile;
+        unlink $zipfile;
+      }
+      if ( ! -f $sdffile ) {
+        $self->debug("$sdffile could not be downloaded.\n");
+        next;
+      }
+      $self->debug("using file $sdffile\n");
+    } # $dy
+  } # $dx
+
+  $self->debug("running Splat...\n");
+  my @args = (
+    '-d', $SRTM_DIR,
+    '-t', "$id.qth", # transmitter site
+    '-o', "$id.ppm", # output file name
+    '-db', $self->max_loss,   # power level range
+    '-L', '30',     # assumed receiver antenna height - XXX should be config
+    '-R', '10',     # maximum coverage radius - XXX should be config
+    '-ngs',         # do not show terrain
+    '-kml',         # generate a KML also
+    '-N',           # suppress noise
+    '-erp', '0',    # transmitter power (only showing relative path loss, so 0)
+  );
+  $self->debug(join(" ", @args) . "\n");
+  my $rv = system($SPLAT_CMD, @args);
+  if ( $rv != 0 ) {
+    $rv = $rv >> 8;
+    die "splat command failed ($rv): $!";
+  }
+
+  # results now in .kml referencing .ppm.
+  # pull the .kml apart to find the coordinates
+  my $kml = XMLin("$id.kml");
+  my $latlon = $kml->{'Folder'}->{'GroundOverlay'}->{'LatLonBox'};
+  $self->box($latlon);
+  my $image = Image::Magick->new;
+  $image->Read("$id.ppm");
+  $self->image($image);
+  # there's a lot we could do with this (GDAL polygonize, etc.) but for
+  # now just hold onto the output.
+}
+
+=item png
+
+Returns the rendered map, as a PNG. The path loss data will be in the red
+channel, scaled according to the I<db_colors> mapping.
+
+=cut
+
+sub png {
+  my $self = shift;
+  my $image = $self->image or return;
+  my $png = $image->Clone;
+  $png->Set(magick => 'png');
+  $png->Transparent('white');
+  return $png->ImageToBlob;
+}
+
+=item mask
+
+Returns the map as a PNG, with solid red in the areas where the path loss is
+below I<max_loss>. The rest of the image will be transparent.
+
+=cut
+
+sub mask {
+  my $self = shift;
+  my $image = $self->image or return;
+  my $png = $image->Clone;
+  $png->Set(magick => 'png');
+  $png->Threshold(threshold => 254, channel => 'Red');
+  $png->Transparent('white');
+  return $png->ImageToBlob;
+}
+
+=item box_geometry
+
+Returns the map area box as a GeoJSON-style geometry object.
+
+=cut
+
+sub box_geometry {
+  my $self = shift;
+  my $box = $self->box;
+  {
+    type => 'Polygon',
+    coordinates => [
+      [ $box->{west}, $box->{south} ],
+      [ $box->{east}, $box->{south} ],
+      [ $box->{east}, $box->{north} ],
+      [ $box->{west}, $box->{north} ],
+      [ $box->{west}, $box->{south} ],
+    ],
+  };
+}
+
+#internal method; you can subclass and override this to create a different
+#antenna pattern.
+
+sub antenna_pattern {
+
+  my $self = shift;
+  my $beamwidth = $self->h_width;
+
+  my $falloff = 2;
+
+  my %az = ();
+  my %el = ();
+  my $theta = 0;
+  while ($theta < 360) {
+
+    if ($theta <= ($beamwidth / 2)) {
+      $az{$theta} = 1.0;
+    } elsif ($theta >= (360 - ($beamwidth / 2))) {
+      $az{$theta} = 1.0;
+    } elsif ($theta > ($beamwidth / 2) and $theta < 180) {
+      my $x = ($theta - ($beamwidth / 2)) / ($beamwidth / 2);
+      $az{$theta} = 1+(-($x*$falloff)**2);
+      $az{$theta} = ($az{$theta} > 0) ? $az{$theta} : 0.001000;
+    } elsif ($theta < (360 - ($beamwidth / 2)) and $theta > 180) {
+      my $x = ($theta - (360 - ($beamwidth / 2))) / ($beamwidth / 2);
+      $az{$theta} = 1+(-(-$x*$falloff)**2);
+      $az{$theta} = ($az{$theta} > 0) ? $az{$theta} : 0.001000;
+    } else {
+      $az{$theta} = 0.001000;
+    }
+
+    $theta++;
+
+  }
+
+  $beamwidth = $self->v_width;
+  $theta = -10;
+  while ($theta <= 90) {
+
+    if (abs($theta) <= ($beamwidth / 2)) {
+      $el{$theta} = 1.0;
+    } else {
+      my $x = (abs($theta) - ($beamwidth / 2)) / ($beamwidth / 2);
+      $el{$theta} = 1+(-($x*$falloff)**2);
+      $el{$theta} = ($el{$theta} > 0) ? $el{$theta} : 0.001000;
+    }
+
+    $theta++;
+  }
+
+#  while ($theta <= 90) {
+#
+#    if ($theta >= -10 and $theta <= 90) {
+#      $el{$theta} = 1.0;
+#    } else {
+#      $el{$theta} = 0.0;
+#    }
+#
+#    $theta += 0.5;
+#
+#  }
+
+  return({%az}, {%el});
+
+}
+
+=back
+
+=head1 AUTHOR
+
+Mark Wells, C<< <mark at freeside.biz> >>
+
+=head1 BUGS
+
+Currently only works in North America due to a hardcoded download path.
+
+The antenna radiation patterns are somewhat arbitrary at the moment.
+
+Please report any bugs or feature requests to C<bug-map-splat at rt.cpan.org>,
+or through the web interface at
+L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Map-Splat>.  I will be
+notified, and then you'll automatically be notified of progress on your bug
+as I make changes.
+
+=head1 SUPPORT
+
+You can find documentation for this module with the perldoc command.
+
+    perldoc Map::Splat
+
+=head1 ACKNOWLEDGEMENTS
+
+John Magliacane (KD2BD) created Splat! and released it under the GNU GPL.
+
+This library is based in part on code by Kristian Hoffmann of Fire2Wire
+Internet Services.
+
+=head1 LICENSE AND COPYRIGHT
+
+Copyright 2016 Mark Wells.
+
+This program is free software; you can redistribute it and/or modify it
+under the terms of the the Artistic License (2.0). You may obtain a
+copy of the full license at:
+
+L<http://www.perlfoundation.org/artistic_license_2_0>
+
+=cut
+
+1;