forked from MapServer/MapServer
-
Notifications
You must be signed in to change notification settings - Fork 2
PerlMapScriptExamples35ex9
Thomas Bonfort edited this page Apr 6, 2012
·
2 revisions
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.
#!perl
#!/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 PerlMapScrip