Splat.pm

   1 package Map::Splat;
   2
   3 use 5.006;
   4 use strict;
   5 use warnings FATAL => 'all';
   6
   7 use File::chdir;
   8 use File::Temp qw(tempdir);
   9 use File::Path qw(make_path);
  10 use File::Slurp qw(read_file write_file);
  11 use XML::LibXML::Simple qw(XMLin);
  12 use List::Util qw(min max);
  13 use Data::Dumper;
  14 use Image::Magick;
  15 use Moo;
  16 use POSIX ('floor', 'ceil');
  17 use v5.10;
  18
  19 =head1 NAME
  20
  21 Map::Splat - Wrapper for the Splat! RF propagation modeling tool
  22
  23 =head1 VERSION
  24
  25 Version 0.01
  26
  27 =cut
  28
  29 our $VERSION = '0.01';
  30
  31 =head1 SYNOPSIS
  32
  33     use Map::Splat;
  34
  35     my $foo = Map::Splat->new();
  36     ...
  37
  38 =head1 METHODS
  39
  40 =over 4
  41
  42 =cut
  43
  44 # requires:
  45 # - Splat
  46 # - wget
  47 # - zip
  48 # - gdal-bin package (not used yet, but just wait)
  49 # - XML::LibXML
  50 # - ImageMagick
  51
  52 my $SRTM_DIR = $ENV{HOME} . '/.splat';
  53 my $SRTM_PATH = 'http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America';
  54 my $SPLAT_CMD = '/usr/bin/splat';
  55
  56 =item new PARAMETERS
  57
  58 Constructs a new mapping scenario. This takes the following parameters:
  59
  60 lon: the tower longitude (east is positive)
  61
  62 lat: the tower latitude (north is positive)
  63
  64 height: the tower height above ground level, in feet
  65
  66 freq: the frequency, in MHz
  67
  68 h_width: the horizontal beamwidth, in degrees
  69
  70 v_width: the vertical beamwidth, in degrees
  71
  72 azimuth: the antenna orientation, in degrees from north (90 = east)
  73
  74 tilt: the downward tilt of the antenna, in degrees
  75
  76 and optionally accepts:
  77
  78 dir: the working directory to use; defaults to a random tempdir
  79
  80 id: the filename prefix for working files; defaults to '1'
  81
  82 name: a string describing the tower; defaults to 'Tower'
  83
  84 min_loss: the lowest path loss level to show on the map; defaults to 80 dB
  85
  86 max_loss: the highest path loss level to show on the map; defaults to 160 dB
  87
  88 =cut
  89
  90 foreach ( qw( lon lat height freq h_width v_width azimuth tilt ) ) {
  91   has $_ => ( is => 'rw', required => 1 );
  92 }
  93
  94 has 'dir' => (
  95   is => 'ro',
  96   default => sub {
  97     tempdir(CLEANUP => 1)
  98   },
  99 );
 100
 101 has 'id' => (
 102   is => 'ro',
 103   default => sub { 1 }
 104 );
 105
 106 has 'name' => (
 107   is => 'ro',
 108   default => sub { 'Tower' }
 109 );
 110
 111 # there are more sophisticated ways to calculate this. for now we just expose
 112 # the ends of the scale as optional parameters.
 113 has 'min_loss' => (
 114   is => 'ro',
 115   default => sub { 80 }
 116 );
 117
 118 has 'max_loss' => (
 119   is => 'ro',
 120   default => sub { 160 }
 121 );
 122
 123 =item image
 124
 125 After calling I<calculate> this will contain an L<Image::Magick> object with
 126 the raw image data.
 127
 128 =cut
 129
 130 has 'image' => (
 131   is => 'rw',
 132 );
 133
 134 =item box
 135
 136 After calling I<calculate> this will contain a hashref of 'north', 'south',
 137 'east', and 'west' indicating the bounds of the map area.
 138
 139 =cut
 140
 141 has 'box' => (
 142   is => 'rw',
 143   default => sub { +{} }
 144 );
 145
 146 =item db_colors
 147
 148 Returns a hashref of (path loss level in dB) => (color level).
 149
 150 =cut
 151
 152 sub db_colors {
 153   my $self = shift;
 154   my %color;
 155   # we have 32 color bands; find the width of a single band
 156   my $db_step = ($self->max_loss - $self->min_loss) / 32;
 157   for (my $i = 0; $i < 32; $i++) {
 158     # + 1 here so that the highest band is _above_ the -db switch setting
 159     my $db = ceil($self->max_loss - ($i * $db_step)) + 1;
 160     my $color = 248 - ($i * 8);
 161     $color{$db} = $color;
 162   }
 163   \%color;
 164 }
 165
 166 sub ew {
 167   my $lon = shift;
 168   $lon < 0 ? 'W' : 'E';
 169 }
 170
 171 sub ns {
 172   my $lat = shift;
 173   $lat < 0 ? 'S' : 'N';
 174 }
 175
 176 sub debug {
 177   my $self = shift;
 178   my $arg = shift;
 179   if (ref $arg) {
 180     warn Dumper($arg);
 181   } else {
 182     warn "$arg\n";
 183   }
 184 }
 185
 186 =item calculate
 187
 188 Invoke Splat to create a map. Currently, the map will cover the 1-degree
 189 quadrangle containing the tower location.
 190
 191 =cut
 192
 193 sub calculate {
 194   my $self = shift;
 195
 196   my $signed_lat = $self->lat;
 197   my $signed_lon = $self->lon;
 198
 199   my $id = $self->id;
 200
 201   local $CWD = $self->dir;
 202
 203   # write the site file
 204   # note that Splat wants longitude in degrees west, from 0 to <360
 205   # (this will also matter for SDF files)
 206   my $splat_lon = -1 * $signed_lon;
 207   $splat_lon += 360 if $splat_lon < 0;
 208
 209   $self->debug("creating QTH file");
 210   write_file("$id.qth",
 211     $self->name . "\n" .
 212     $signed_lat . "\n" .
 213     $splat_lon . "\n" .
 214     $self->height . "\n"
 215   );
 216
 217   # Longley-Rice parameters
 218   my $freq = $self->freq;
 219   my $pol = 1; # polarization = vertical
 220   $self->debug("creating LRP file for $freq MHz");
 221
 222   my $lrp = <<EOF;
 223 15.000  ; Earth Dielectric Constant (Relative permittivity)
 224 0.005   ; Earth Conductivity (Siemens per meter)
 225 301.000 ; Atmospheric Bending Constant (N-Units)
 226 $freq   ; Frequency in MHz (20 MHz to 20 GHz)
 227 5       ; Radio Climate
 228 $pol    ; Polarization (0 = Horizontal, 1 = Vertical)
 229 0.50    ; Fraction of situations
 230 0.50    ; Fraction of time
 231 0   ; ERP
 232 EOF
 233   write_file("$id.lrp", $lrp);
 234
 235   # Loss contour file: specifies the color mapping for the output
 236   $self->debug("creating $id.lcf\n");
 237   my $lcf = '';
 238   my $colors = $self->db_colors;
 239   for my $db (sort { $a <=> $b } keys %$colors) {
 240     my $red = $colors->{$db};
 241     $lcf .= "$db: $red, 0, 0\n";
 242   }
 243   write_file("$id.lcf", $lcf);
 244
 245   # AZ file: radiation pattern in horizontal plane
 246   # should be able to take this as an input (and kinda can, by overriding
 247   # antenna_pattern)
 248   my $bearing = $self->azimuth;
 249   my $downtilt = $self->tilt;
 250   my ($az, $el) = $self->antenna_pattern;
 251   $self->debug("creating $id.az for bearing $bearing\n");
 252
 253   open my $fh, '>', "$id.az";
 254   printf $fh "%0.1f\n", ($bearing);
 255   foreach my $t (sort { $a <=> $b } keys (%$az)) {
 256     printf $fh "%d\t%0.6f\n", ($t, $az->{$t});
 257   }
 258
 259   close $fh;
 260
 261   # EL file: radiation pattern in vertical plane
 262   $self->debug("creating $id.el for downtilt $downtilt\n");
 263   open $fh, '>', "$id.el";
 264   printf $fh "%0.1f\t%0.1f\n", ($downtilt, $bearing);
 265   foreach my $t (sort { $a <=> $b } keys (%$el)) {
 266     printf $fh "%0.1f\t%0.3f\n", ($t, $el->{$t});
 267   }
 268   close $fh;
 269
 270   # make sure we have SRTM for the area around the antenna
 271   # ("area around the antenna" should be configurable)
 272   # (and support SRTM1/high def version? only works in North America)
 273
 274   for (my $dx = -2; $dx <= 2; $dx++) {
 275     for (my $dy = -2; $dy <= 2; $dy++) {
 276       my $splat_lon_east = floor($splat_lon) + $dx;
 277       my $splat_lon_west = ($splat_lon_east + 1) % 360;
 278       my $lat_south = floor($signed_lat) + $dy;
 279       my $lat_north = $lat_south + 1;
 280
 281       next if $lat_north > 90 or $lat_south < -90; # haha
 282
 283       # the naming convention used by Splat's SDF files:
 284       # south:north:east:west
 285       # longitude is degrees west, from 0 to 359.
 286       # note that the last one (from 1E to 0W) is 359:0, not 359:360
 287       my $sdffile = "$SRTM_DIR/$lat_south:$lat_north:$splat_lon_east:$splat_lon_west.sdf";
 288       if ( ! -f $sdffile ) {
 289         if ( ! -d $SRTM_DIR ) {
 290           make_path($SRTM_DIR); # dies on error
 291         }
 292         local $CWD = $SRTM_DIR;
 293         # the USGS, though, just marks them as N/S and E/W. reference is the
 294         # southwest corner.
 295         # that's the negative of splat longitude, shifted to the range (-180, 180)
 296         # or in the east, it's 360 - splat longitude
 297         my $usgs_lon_west = -1 * $splat_lon_west;
 298         $usgs_lon_west += 360 if $usgs_lon_west < -180;
 299         my $hgtfile = ns($lat_south) . abs($lat_south) .
 300                       ew($usgs_lon_west) . abs($usgs_lon_west) .
 301                       '.hgt';
 302         my $zipfile = $hgtfile . '.zip';
 303         $self->debug("fetching $SRTM_PATH/$zipfile\n");
 304         my $status = system('wget', "$SRTM_PATH/$zipfile");
 305         if ( $status > 0 ) {
 306           # normal if the tower location is near a coast, so try to continue
 307           warn "couldn't download $zipfile\n";
 308           next;
 309         }
 310         $status = system('unzip', $zipfile);
 311         # this isn't normal though
 312         die "couldn't extract $zipfile\n" if $status > 0;
 313         # unfortunately gives no useful output
 314         $status = system('srtm2sdf', $hgtfile);
 315         unlink $hgtfile;
 316         unlink $zipfile;
 317       }
 318       if ( ! -f $sdffile ) {
 319         $self->debug("$sdffile could not be downloaded.\n");
 320         next;
 321       }
 322       $self->debug("using file $sdffile\n");
 323     } # $dy
 324   } # $dx
 325
 326   $self->debug("running Splat...\n");
 327   my @args = (
 328     '-d', $SRTM_DIR,
 329     '-t', "$id.qth", # transmitter site
 330     '-o', "$id.ppm", # output file name
 331     '-db', $self->max_loss,   # power level range
 332     '-L', '30',     # assumed receiver antenna height - XXX should be config
 333     '-R', '10',     # maximum coverage radius - XXX should be config
 334     '-ngs',         # do not show terrain
 335     '-kml',         # generate a KML also
 336     '-N',           # suppress noise
 337     '-erp', '0',    # transmitter power (only showing relative path loss, so 0)
 338   );
 339   $self->debug(join(" ", @args) . "\n");
 340   my $rv = system($SPLAT_CMD, @args);
 341   if ( $rv != 0 ) {
 342     $rv = $rv >> 8;
 343     die "splat command failed ($rv): $!";
 344   }
 345
 346   # results now in .kml referencing .ppm.
 347   # pull the .kml apart to find the coordinates
 348   my $kml = XMLin("$id.kml");
 349   my $latlon = $kml->{'Folder'}->{'GroundOverlay'}->{'LatLonBox'};
 350   $self->box($latlon);
 351   my $image = Image::Magick->new;
 352   $image->Read("$id.ppm");
 353   $self->image($image);
 354   # there's a lot we could do with this (GDAL polygonize, etc.) but for
 355   # now just hold onto the output.
 356 }
 357
 358 =item png
 359
 360 Returns the rendered map, as a PNG. The path loss data will be in the red
 361 channel, scaled according to the I<db_colors> mapping.
 362
 363 =cut
 364
 365 sub png {
 366   my $self = shift;
 367   my $image = $self->image or return;
 368   my $png = $image->Clone;
 369   $png->Set(magick => 'png');
 370   $png->Transparent('white');
 371   return $png->ImageToBlob;
 372 }
 373
 374 =item mask
 375
 376 Returns the map as a PNG, with solid red in the areas where the path loss is
 377 below I<max_loss>. The rest of the image will be transparent.
 378
 379 =cut
 380
 381 sub mask {
 382   my $self = shift;
 383   my $image = $self->image or return;
 384   my $png = $image->Clone;
 385   $png->Set(magick => 'png');
 386   $png->Threshold(threshold => 254, channel => 'Red');
 387   $png->Transparent('white');
 388   return $png->ImageToBlob;
 389 }
 390
 391 =item box_geometry
 392
 393 Returns the map area box as a GeoJSON-style geometry object.
 394
 395 =cut
 396
 397 sub box_geometry {
 398   my $self = shift;
 399   my $box = $self->box;
 400   {
 401     type => 'Polygon',
 402     coordinates => [
 403       [ $box->{west}, $box->{south} ],
 404       [ $box->{east}, $box->{south} ],
 405       [ $box->{east}, $box->{north} ],
 406       [ $box->{west}, $box->{north} ],
 407       [ $box->{west}, $box->{south} ],
 408     ],
 409   };
 410 }
 411
 412 #internal method; you can subclass and override this to create a different
 413 #antenna pattern.
 414
 415 sub antenna_pattern {
 416
 417   my $self = shift;
 418   my $beamwidth = $self->h_width;
 419
 420   my $falloff = 2;
 421
 422   my %az = ();
 423   my %el = ();
 424   my $theta = 0;
 425   while ($theta < 360) {
 426
 427     if ($theta <= ($beamwidth / 2)) {
 428       $az{$theta} = 1.0;
 429     } elsif ($theta >= (360 - ($beamwidth / 2))) {
 430       $az{$theta} = 1.0;
 431     } elsif ($theta > ($beamwidth / 2) and $theta < 180) {
 432       my $x = ($theta - ($beamwidth / 2)) / ($beamwidth / 2);
 433       $az{$theta} = 1+(-($x*$falloff)**2);
 434       $az{$theta} = ($az{$theta} > 0) ? $az{$theta} : 0.001000;
 435     } elsif ($theta < (360 - ($beamwidth / 2)) and $theta > 180) {
 436       my $x = ($theta - (360 - ($beamwidth / 2))) / ($beamwidth / 2);
 437       $az{$theta} = 1+(-(-$x*$falloff)**2);
 438       $az{$theta} = ($az{$theta} > 0) ? $az{$theta} : 0.001000;
 439     } else {
 440       $az{$theta} = 0.001000;
 441     }
 442
 443     $theta++;
 444
 445   }
 446
 447   $beamwidth = $self->v_width;
 448   $theta = -10;
 449   while ($theta <= 90) {
 450
 451     if (abs($theta) <= ($beamwidth / 2)) {
 452       $el{$theta} = 1.0;
 453     } else {
 454       my $x = (abs($theta) - ($beamwidth / 2)) / ($beamwidth / 2);
 455       $el{$theta} = 1+(-($x*$falloff)**2);
 456       $el{$theta} = ($el{$theta} > 0) ? $el{$theta} : 0.001000;
 457     }
 458
 459     $theta++;
 460   }
 461
 462 #  while ($theta <= 90) {
 463 #
 464 #    if ($theta >= -10 and $theta <= 90) {
 465 #      $el{$theta} = 1.0;
 466 #    } else {
 467 #      $el{$theta} = 0.0;
 468 #    }
 469 #
 470 #    $theta += 0.5;
 471 #
 472 #  }
 473
 474   return({%az}, {%el});
 475
 476 }
 477
 478 =back
 479
 480 =head1 AUTHOR
 481
 482 Mark Wells, C<< <mark at freeside.biz> >>
 483
 484 =head1 BUGS
 485
 486 Currently only works in North America due to a hardcoded download path.
 487
 488 The antenna radiation patterns are somewhat arbitrary at the moment.
 489
 490 Please report any bugs or feature requests to C<bug-map-splat at rt.cpan.org>,
 491 or through the web interface at
 492 L<http://rt.cpan.org/NoAuth/ReportBug.html?Queue=Map-Splat>.  I will be
 493 notified, and then you'll automatically be notified of progress on your bug
 494 as I make changes.
 495
 496 =head1 SUPPORT
 497
 498 You can find documentation for this module with the perldoc command.
 499
 500     perldoc Map::Splat
 501
 502 =head1 ACKNOWLEDGEMENTS
 503
 504 John Magliacane (KD2BD) created Splat! and released it under the GNU GPL.
 505
 506 This library is based in part on code by Kristian Hoffmann of Fire2Wire
 507 Internet Services.
 508
 509 =head1 LICENSE AND COPYRIGHT
 510
 511 Copyright 2016 Mark Wells.
 512
 513 This program is free software; you can redistribute it and/or modify it
 514 under the terms of the the Artistic License (2.0). You may obtain a
 515 copy of the full license at:
 516
 517 L<http://www.perlfoundation.org/artistic_license_2_0>
 518
 519 =cut
 520
 521 1;