tcheck.pl

The cl.tar.gz url is ​http://www.highwayengineer.co.medina.oh.us/cl.tar.gz , unfortunately it is large for modem download but hopefully it will be useful to other users.

#!/usr/bin/perl
#
# Copyright (C) 2002, Lowell Filak.
# You may distribute this file under the terms of the Artistic
# License.
#
# Given the name of an existing point shapefile, the name of an existing
#   line shapefile, a radius distance, a field name in the point shapefile,
#   a filed name in the line shapefile, & optionally an accuracy field name
#   in the point shapefile this routine will search within a circular buffer
#   of each point (of specified radius + optional accuracy) for a line with a matching
#   field value and then assign the point a 0 error for a single match. Otherwise the point
#   will be assigned an error of 1 for no matches & x for a number of
#   matches > 1 where x is the number of matches.
#   Note: Currently the field matching portion contains syntax to reduce
#         the value of the fields to numeric before matching.

#         The routine only works if both shapefiles are in the current dir.
#         The routine requires a RECORD field in the point & line shapefiles that
#           reflects the number of the record (this must line up) with the
#           shape number of the point. First record number being 0.
#         Instead of using the xbase module directly like find.pl this
#           routine uses the sql style interface to the shapefile attributes.
#           Why? Because it is a lot easier to code.
#         The routine requires a field named errflag to exist in the point
#           shapefile (dbf).
#         The errflag field created by tcounts.pl is limited to a max of 99.
#           This means that the point buffer should be small enough not to
#           overlap more than 99 lines.
#
# Required modules are mapscript (installed as part of make install
#   http://mapserver.gis.umn.edu),
#   Getopt (normally included with Perl),
#   DBI (cpan), DBD::XBase (cpan),
#   & XBase (cpan).
#   Please run tcounts.pl first then:
#     please download cl.tar.gz also, and:
#     tar -xf cl.tar.gz --ungzip
#
# Suggested run line = ./tcheck.pl -pfile=traffic -lfile=cl -iradius=35 -pfield=roadnumber -lfield=first_ -afield=accuracyt
#
# Include the mapscript module.
use mapscript;
#
# Include the xbase and dbi modules for searching and updating values.
use XBase;
use DBI;
#
# Include the getopt module to read input.
use Getopt::Long;
#
# Grab the file name from the input.
&GetOptions('pfile=s' => \$pfile, 'lfile=s' => \$lfile, 'iradius=i' => \$iradius, 'pfield=s' => \$pfield, 'lfield=s' => \$lfield, 'afield=s' => $afield);
if ( !$pfile || !$lfile || !$iradius || !$pfield || !$lfield) {
  print "Syntax: tcheck.pl -pfile=[in_point_shapefile_name] -lfile=[in_line_shapefile_name] -iradius=[search_radius] -pfield=[point_shapefile_field_name] -lfield=[line_shapefiel_field_name] -afield={optional_accuracy_field_name_in_point_shapefile";
  exit 0;
}
#
# Create a unique name for a new mapfile for querying the road centerlines.
#
# Grab the date.
my ($sec,$min,$hr,$mdy,$mnth,$yr,$wdy,$ydy,$isdst) = localtime;
#
# Grab the process id.
my $spid = $$;
#
# Create the name & make sure it is no longer than 8 characters.
my $sfile = "T$hr$min$sec$spid";
#my $sfile = "TEST";
#
# Open the mapfile for writing.
open(MAPFILE, ">$sfile.map");
#
# Open the existing point shapefile.
my $inshpf = new shapefileObj($pfile, -1) or die "Unable to open shapefile $pfile.";
#
# What are the extents.
my $inshpminx = $inshpf->{bounds}->{minx};
my $inshpminy = $inshpf->{bounds}->{miny};
my $inshpmaxx = $inshpf->{bounds}->{maxx};
my $inshpmaxy = $inshpf->{bounds}->{maxy};
#
# Create the contents of the mapfile.
print MAPFILE <<EOF;
#
NAME $sfile
STATUS ON
SIZE 600 600
SYMBOLSET "$sfile.sym"
EXTENT $inshpminx $inshpminy $inshpmaxx $inshpmaxy
UNITS FEET
SHAPEPATH ""
IMAGECOLOR 255 255 255
LAYER
  NAME centerline
  TYPE LINE
  STATUS ON
  DATA "$lfile"
  TEMPLATE 'bogus.html'
  CLASS
    COLOR 255 0 0
    NAME "Centerline"
  END
END 
END
EOF
#
# Close the mapfile.
close MAPFILE;
#
# Open the symbol file for writing.
open(SYMFILE, ">$sfile.sym");
#
# Create the contents of the symbol file.
print SYMFILE <<EOF;
SYMBOLSET
END
EOF
#
# Close the symbol file.
close SYMFILE;
#
# How many points are there.
my $innumshp = $inshpf->{numshapes};
#
# Connect to the dbf files as if they were an rdbms.
$dbf=DBI->connect("dbi:XBase:.");
#
# Create a circle shape.
#   Thanks to Tim Mackey.

#
# Create a point object to hold the line-segment increments.
my $p=new pointObj();
#
# Set the value for pi.
my $pi=3.141592654;
#
# Create the mapfileobject.
my $map = new mapObj("$sfile.map") or die("Unable to Open Default MapFile $sfile.map!");
#
# Create the layer object for the later queries.
my $lyr = $map->getLayerByName("centerline") or die('Unable to Open Centerline Layer!');
#

# Set the default query result to blank.
my @row = ();
#
# Set the starting radius to 0.
my $radius = 0;
#
# Create the point holder.
my $inpnt = new pointObj();
#
# Loop through each point.
for ($inpntnum=0; $inpntnum<$innumshp; $inpntnum++ ) {
  print "Checking Point #$inpntnum...\n";
  #
  # Grab the point by number.
  my $junk = $inshpf->getPoint($inpntnum, $inpnt);
  #
  # Set the default accuracy to 0.
  my $accuracy = 0;

  #
  # If the accuracy field is specified use it.
  if ( $afield ) {
    #
    # Retrieve the accuracy distance for this point by point number.
    $sth = $dbf->prepare("SELECT $afield FROM $pfile WHERE record = $inpntnum");
    $sth->execute;
    #
    # Grab the result of the query.
    @row = $sth->fetchrow_array;
    #
    # Set the value for the accuracy.
    $accuracy = $row[0];
  }
  #
  # Set the value for the radius.
  $radius = $accuracy + $iradius;
  #
  # What is the field value for this point.
  $sth = $dbf->prepare("SELECT $pfield FROM $pfile WHERE record = $inpntnum");
  $sth->execute;

  #
  # Grab the result of the query.
  @row = $sth->fetchrow_array;
  #
  # Set the value for the field.
  $pfldval = $row[0];
  #
  # Set the coordinates of the center of the circle to the coordinates of the
  #   existing point.
  my $x = $inpnt->{x};
  my $y = $inpnt->{y};
  #
  # Create a line object to hold the bounds of the shape.
  my $line=new lineObj();
  #
  # Create the circle shapeobject.
  my $circle=new shapeObj($mapscript::MS_SHAPE_POLYGON);
  #
  # Loop through the segments of the circle.
  for($i=0;$i<2.1*$pi;$i=$i+$pi/100) {
    #
    # The x & y coordinates of the segment point.
    $p->{x} = $x + $radius * cos($i);
    $p->{y} = $y + $radius * sin($i);
    #
    # Add the point to the line.
    $line->add($p);
    #
    # Set a marker for the first/last point if this is the first point.
    # This is not needed. See note below.
    #if ( $i==0 ) {
      #$firstpntx = $p->{x};
      #$firstpnty = $p->{y};
    #}
  }
  #
  # Throw in the closing point.
  # This must not be needed. The original code never added this
  #   point to the line but yet the shape closes.?
  #$p->{x} = $firstpntx;
  #$p->{y} = $firstpnty;
  #
  # Add the line to form the circle shape.
  $circle->add($line);
  #
  # Clear out the line object.
  undef $line;
  #
  # Query the line layer with this shape.
  $lyr->queryByShape($map, $circle);
  #
  # Clear out the circle object.
  undef $circle;
  #
  # How many matches found.
  #
  # Create a resultcache object to see how many results.
  my $rsltcache = $lyr->{resultcache};
  #
  # How many matches did we find.
  my $numrslt = $rsltcache->{numresults};
  #
  # Clear the resultcache.
  undef $rsltcache;
  #
  # if there is no match then mark point with error of 1.
  if ( $numrslt == 0 ) {
    $sth = $dbf->prepare("UPDATE $pfile SET errflag = 1 WHERE record = $inpntnum");
    $sth->execute;
  }
  #
  # Set the match quantity to the number of possible matches to start with.
  my $possmatches = $numrslt;
  #
  # Loop through each match to see if the field value matches the point.
  for ($line=0; $line<$numrslt; $line++) {
    #
    # Grab the nth result of the query.
    $resultmember = $lyr->getResult($line);    
    #
    # What is the shape number of this possible match.
    $shaperecnum = $resultmember->{shapeindex};
    #
    # Query the line attribute to see if the fields match.
    $sth = $dbf->prepare("SELECT $lfield FROM $lfile WHERE record = $shaperecnum");
    $sth->execute;
    #
    # Grab the result of the query.

    @row = $sth->fetchrow_array;
    #
    # Set the value of the field.
    $lfldval = $row[0];
    #
    # Take out anything that isn't an integer.
    $lfldval =~ s/[^\x30-\x39]//g;
    #
    # Does this line value match the point value.
    if ( $pfldval == $lfldval ) {
      #
      # Yes it does.
      #
      # I propose to do nothing.
    }
     else {
      #
      # Nope.
      #
      # Subtract one from the possible matches.
      $possmatches = $possmatches - 1;
     }
  }
  #
  # If there is no match found then mark point with error.
  if ( $possmatches == 0 ) {
    #
    # Mark with a 1.
    print "\tNO Match Found...\n";
    $sth = $dbf->prepare("UPDATE $pfile SET errflag = 1 WHERE record = $inpntnum");
    $sth->execute;
  }
  elsif ( $possmatches > 1 ) {
    #
    # Mark with the number of matches.
    print "\t$possmatches Matches Found...\n";
    $sth = $dbf->prepare("UPDATE $pfile SET errflag = $possmatches WHERE record = $inpntnum");
    $sth->execute;
  }
  elsif ( $possmatches == 1 ) {
    #
    # Make sure the error is set to 0 just incase the routine was run multiple
    #   times on the same shapefile with different radii.
    print "\tOne Match Found...\n";
    $sth = $dbf->prepare("UPDATE $pfile SET errflag = 0 WHERE record = $inpntnum");
    $sth->execute;
  }
}   
#
# Get rid of the temporary mapfile & symbol file.
unlink "$sfile.map";

unlink "$sfile.sym";

---

back to PerlMapScript