path: root/Splat.pm

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;