wiki:PerlMapScriptExamples35ex8

tcounts.pl

The url for the data is http://www.highwayengineer.co.medina.oh.us/StickeDB.pdb .

#!/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 palm database and the name
#   of an existing point shapefile this routine will append the points
#   in the pdb to the shapefile.
# Given just the name of an existing shapefile the routine will attempt
#   to use pilot-xfer to download the pdb file and append the points to the
#   shapefile.
# Given just the name of an existing pdb file this routine will create a
#   new shapefile of the points in the pdb file.
# Given neither name this routine will attempt to use pilot-xfer to download
#   the pdb file and create a new point shapefile from it.
# Notes: The fields in the pdb file should match the fields in the
#        existing shapefile (dbf) or the assignments will either be
#        wrong or will cause the routine to bomb.
#   The pdb file for download is assumed to be StickeDB.pdb and the routine
#     is written to read sticke style pdb files only.
#   The pdb reading section of this routine is not complete but is setup
#     somewhat generic and should be extendable to any sticke database
#     schema.
#   The pilot-xfer download line assumes the default pilot device
#     (/dev/pilot) exists.
#   The routine also assumes a *n*x system, please change command
#     lines accordingly.
#   If nad support is needed in proj.4 please verify that the ntv1_can.dat
#     file is included before compiling. If not, grab a newer release.
#
# Required modules are mapscript (installed as part of make install
#   http://mapserver.gis.umn.edu),
#   Getopt (normally included with Perl),
#   Palm (p5-Palm-1.1.5 http://theoryx5.uwinnipeg.ca/CPAN/data/p5-Palm/Palm/Raw.html),
#   & XBase (cpan).
#   Please download StickeDB.pdb also.
#
# Additional requirements are: a working pilot-xfer (pilot-link http://www.pilot-link.org/)
#   installation,
#   a working StickePad.prc and StickePlates.prc (StickeV2Programs
#   http://www.cs.ukc.ac.uk/projects/mobicomp/Fieldwork/) on a
#   PalmOS handheld device, & a working proj.4 install compiled with
#   the optional nad files in place and with the cs2cs
#   command working (www.remotesensing.org/proj/).
#
# Current GPS information: Palm IIIX-PDA, Garmin Etrex Summit-GPS, Blue
#   Hills Innovations-Garmin2Palm cable (http://www.blue-hills-innovations.com).
#
# Suggested run line = ./tcounts.pl -pdbfile=StickeDB -sfile=traffic
#
# Syntax: tcounts.pl -pdbfile=[in_pdb_filename] -sfile=[out_shapefile_name]
#
# Include the pdb and pdb-raw modules.
use Palm::PDB;
use Palm::Raw;
#
# Include the mapscript module.
use mapscript;
#
# Include the xbase module for creating the dbf records.
use XBase;
#
# Include the getopt module to read input.
use Getopt::Long;
#
# Helpful definitions for StickeDB.pdb:
#   I view the structure as very similar to an rdbms.
#   Database - refers to the pdb file itself.
#   Table - refers to the rdbms-like table name included on every record.
#     Note: Each record can belong to a different table or even a different
#           table deffinition under the same table name!
#     Note: Tables are refered to as 'templates' inside sticke.
#   Record - refers to the entire 'data' portion returned by the pdb->data obj.
#     Note: Records are refered to as 'notes' inside of sticke.
#   Field - refers to the section of the 'data' portion of the record which
#           spans from the beginning of one field name to the beginning of the
#           next.
#     Note: Fields are refered to as both 'fields' and 'items' inside of sticke.
#   Part - refers to the sections of a field that define the schema of the
#          field (schema - data type, constraints, etc).
#
# Grab the file names from the input.
&GetOptions('pdbfile=s' => \$pdbfile, 'sfile=s' => \$sfile);
if ( !$sfile ) {
  #
  # Create a unique name for the shapefile based on date and process number.

  #
  # Grab the time.
  ($sec,$min,$hr,$mdy,$mnth,$yr,$wdy,$ydy,$isdst) = localtime;
  #
  # Grab the process id.
  $spid = $$;

  #
  # Create the name & make sure it is no longer than 8 characters.
  $sfile = substr("$hr$min$sec$spid", -8);
}
if( !$pdbfile ) {
  #
  # Download the pdb file.
  system("pilot-xfer -f StickeDB");
  #
  # Set the pdb file name.
  $pdbfile = 'StickeDB';
}
#
# Create the pdb object on the pdb file.
my $pdb = new Palm::PDB;
$pdb->Load("$pdbfile.pdb");
#
# How many pdb records are there.
my @records = @{$pdb->{records}};
my $numrecs = scalar(@records);
#
# Create the information array for each field/data type.
#   Note: The routine is incomplete at this time because we currently only use
#    7-number & 4-location but all types are here for documentation.
my @types = ('bearing', 'boolean', 'date', 'textline', 'location', 'unused', 'notepad', 'number', 'picklist', 'subnote', 'time');
#
# Create the field type number-of-parts array.
my @parts = (0,0,0,0,23,0,0,11,0,0,0);
#
# Create the field types unpack string array.
#   For location [offset+10]=dontno5, [offset+11]=latdegrees,
#    [offset+12]=dontno6, [offset+13]=latminuteswhole,
#    [offset+14]=dontno7, [offset+15]=latminutesdecimal,
#    [offset+16]=elevation, [offset+17]=dontno8, [offset+18]=longdegrees,
#    [offset+19]=dontno9, [offset+20]=longminuteswhole,
#    [offset+21]=dontno10, [offset+22]=longminutesdecimal
#    Note: At this point I can't find the NSEW/+- indication for lat/lon!
#  For number [offset+10]=number
my @ustring = (' a*', ' a*', ' a*', ' a*', ' l n a2 n a2 n s a2 n a2 n a2 n', ' a*', ' a*', ' N', ' a*', ' a*', ' a*');
#
# Create the array for the pdb-field-type to dbf-field-type conversion.
my @dbfftype = ('C', 'L', 'D', 'C', 'C', 'C', 'C', 'N', 'C', 'C', 'N');
#

# Create the array for the pdb-field-size to dbf-field-size conversion.
my @dbffsize = ('255', '1', '8', '255', '31', '0', '255', '11', '255', '255', '10');
#
# Create the array for the pdb-field-decimal to dbf-field-decimal conversion.
# No need for this now. All decimals should be null going into the dbf file.
#
# Create the array for tracking field offsets.
#   To be used later while creating dbf records for each gps point.
my @offsets = ();
#
# Create the initial unpack string.
my $unpackstr = "";
#
# Create the initial data array.
my @recordinfo = ();
#
# Initialize the dbf record count to 0.
my $dbfreccnt = 0;
#
# Does the dbf file already exist or does it yet to exist.
if ( -e "$sfile.dbf" ) {
  #
  # Open the existing dbf file for appending to.
  $dbh = new XBase "$sfile.dbf" or die XBase->errstr;
  #
  # To be able to increment the record # starting at the last existing
  #   record how many records are there.
  $dbfreccnt = $dbh->last_record + 1;
}
 else {
  $dbfreccnt = -1;
 }

#
# Does the shapefile already exist or is it yet to exist.
if ( -e "$sfile.shp" ) {
  #
  # Move the existing shapefile to a temp name.
  #   This is done because the -1 option on shapefileObj open
  #     only allows for read not append.
  #   However as (hopefully) shown below this is not hard to
  #     implement inside mapscript.
  system("mv $sfile.shp thistemp.shp; mv $sfile.shx thistemp.shx; touch thistemp.dbf");  
  #
  # Open the existing file.
  $ecounts = new shapefileObj("thistemp",-1);
  #
  # Create the replacement shapefile.
  $tcounts = new shapefileObj("$sfile",$mapscript::MS_SHAPEFILE_POINT);
  #
  # Create the transfer point object.
  my $trnspnt = new pointObj();  
  #
  # Loop through each existing point and recreate it in the new shapefile.
  for ($epnt=0; $epnt<$dbfreccnt; $epnt++) {
    #
    # Get the existing point.
    $ecounts->getPoint($epnt,$trnspnt);
    #
    # Put the point into the new shapefile.
    $tcounts->addPoint($trnspnt);
  }
}
 else {
  #
  # Create the new file.
  $tcounts = new shapefileObj("$sfile",$mapscript::MS_SHAPEFILE_POINT);
 }
#
# Create the point object for insertion into the shapefile.
my $pnt = new pointObj();
#
# Loop through each record.
#   [0]=dontno1, [1]=tablenamechars, [2]=tablename, [3]=dontno2, [4]=#offields
for ($recrd=1; $recrd<$numrecs; $recrd++) {
  #
  # Create the array for tracking field offsets.
  #   To be used later while creating dbf records for each gps point.
  my @offsets = ();
  #
  # Set the initial value for unpacking table name.
  $unpackstr = "a38 n";
  @recordinfo = unpack($unpackstr, $records[$recrd]->{data});
  #
  # The character count returned is actual + 1.
  $recordinfo[1] = $recordinfo[1] - 1;
  #
  # Unpack the dontno1 and the table name length.
  @recordinfo = unpack($unpackstr, $records[$recrd]->{data});
  #
  # The character count returned is actual + 1.
  $recordinfo[1] = $recordinfo[1] - 1;

  #
  # Add the remainder of the record info to the unpack string
  #   (name, dontno2, #offields).
  $unpackstr = $unpackstr . " a$recordinfo[1] a3 n";
  #
  # Add the first 10 parts of the first field info to the unpack string.
  #   All fields appear to have these in common even if they are blank
  #    for fields that do not use the particular part.
  #   (fieldname, dontno3, datatype, isrange, null, upperlimit, lowerlimit,
  #    step, dontno4, fieldsize).
  #   [offset+0]=fieldname, [offset+1]=dontno3, [offset+2]=datatype
  #   [offset+3]=isrange, [offset+4]=null,
  #   [offset+5]=uppperlimit, [offset+6]=lowerlimit, [offset+7]=step,
  #   [offset+8]=dontno4, [offset+9]=fieldsize
  $unpackstr = $unpackstr . " A19 a10 n n a N N N a14 n";
  #
  # Set the inital value for the field offset
  #   (the number of parts for the previous field(s)).
  my $fieldoffset = 0;
  #
  # Grab the rest of the record info and
  #   the field info up to the data length indicator.
  ($recordinfo[0], $recordinfo[1], $recordinfo[2], $recordinfo[3], $recordinfo[4], @fieldinfo) = unpack $unpackstr, $records[$recrd]->{data};
  #
  # The character count returned is actual + 1.
  $recordinfo[1] = $recordinfo[1] - 1;
  #
  #
  # Print the record info to see if we got this right.
  #print "\nRecord # = $recrd\nNumber of Characters in Table Name = $recordinfo[1]\nTable Name = $recordinfo[2]\nNumber of Fields = $recordinfo[4]\n";
  # Loop through each field.
  for ($fld=0; $fld<$recordinfo[4]; $fld++) {
    #
    # The actual field number to print is fld + 1.
    my $fldprint = $fld + 1;
    #
    # Grab the field info up to the data length indicator.
    ($recordinfo[0], $recordinfo[1], $recordinfo[2], $recordinfo[3], $recordinfo[4], @fieldinfo) = unpack $unpackstr, $records[$recrd]->{data};
    #
    # The character count returned is actual + 1.
    $recordinfo[1] = $recordinfo[1] - 1;
    #
    # What is the length of the data.

    my $dlength = $fieldinfo[$fieldoffset+9];
    #
    # The field type comes in strange sometimes so this should truncate it
    #   so it only contains values of 0-10.
    $fieldinfo[$fieldoffset+2] = 256 * ( ( $fieldinfo[$fieldoffset+2] / 256 ) - ( int( $fieldinfo[$fieldoffset+2] / 256 ) ) );
    #
    # Okay, the same thing happens with the range.
    $fieldinfo[$fieldoffset+3] = 256 * ( ( $fieldinfo[$fieldoffset+3] / 256 ) - ( int( $fieldinfo[$fieldoffset+3] / 256 ) ) );
    #
    # Add to the unpack string the unpack string for the field type.
    $unpackstr = $unpackstr . $ustring[$fieldinfo[$fieldoffset+2]];
    #
    # For some reason the type appears to be 8-bit instead if 16. So
    #   to make sure 
    #
    # Add to the array the rest of the parts for the field.
    ($recordinfo[0], $recordinfo[1], $recordinfo[2], $recordinfo[3], $recordinfo[4], @fieldinfo) = unpack $unpackstr, $records[$recrd]->{data};
    #
    # Escape out any unprintable characters in the field name.
    $fieldinfo[$offsets[$iname]+0] = uc($fieldinfo[$offsets[$iname]+0]);
    $fieldinfo[$offsets[$iname]+0] =~ s/[^\x41-\x5A]//g;
    #
    # The field type comes in strange sometimes so this should truncate it
    #   so it only contains values of 0-10.
    #   (binary/unpack guru applications now being accepted).
    # Note: Basically this divides by base 16 to move the number 2 decimal
    #       places left then truncates the whole number and multiplies by
    #       base 16 to move the decimal 2 places right.
    $fieldinfo[$fieldoffset+2] = 256 * ( ( $fieldinfo[$fieldoffset+2] / 256 ) - ( int( $fieldinfo[$fieldoffset+2] / 256 ) ) );
    #
    # Okay, the same thing happens with the range.
    $fieldinfo[$fieldoffset+3] = 256 * ( ( $fieldinfo[$fieldoffset+3] / 256 ) - ( int( $fieldinfo[$fieldoffset+3] / 256 ) ) );
    #
    # Print the field info to see if we got this right.
    #print "Field Offset = $fieldoffset\nField $fldprint Name = $fieldinfo[$fieldoffset+0]\nData Type = $fieldinfo[$fieldoffset+2]\nIsRange = $fieldinfo[$fieldoffset+3]\nUpper Limit = $fieldinfo[$fieldoffset+5]\nLower Limit = $fieldinfo[$fieldoffset+6]\nStep = $fieldinfo[$fieldoffset+7]\nField Size = $fieldinfo[$fieldoffset+9]\n";
    #
    # How many data parts are there.
    #   The total number of field parts - 10 is the number of data parts.
    my $dparts = $parts[$fieldinfo[$fieldoffset+2]];
    #
    # Loop through each of the field value parts.
    for ($dpart=10; $dpart<$dparts; $dpart++) {
      #
      # The actual data part id is the current dpart - 9 (0 thru 9 of the
      #   field array).
      my $dprint = $dpart - 9;
      #
      # Print the field info to see if we got this right.
      #print "Data Value $dprint = $fieldinfo[$fieldoffset+$dpart]\n";
    }
    #
    # If the field is a location convert the lat/long to state plane.
    if (  $fieldinfo[$fieldoffset+2] == 4 ) {
      #
      # Do the convert.
      # Bunches of notes: The projection name is latlong but supply
      #   the coordinates as long/lat.
      #   The +to section contains units of us-ft but MUST specify
      #     false_east(x_0) in meters.
      #   An indespensible resource was:
      #     http://www.edc.uri.edu/nrs/classes/NRS522/Tools/StatePlaneZones.htm
      # Note: If I was smart I would have used the pointObj project method.
      system("echo \'$fieldinfo[$fieldoffset+18]d$fieldinfo[$fieldoffset+20].$fieldinfo[$fieldoffset+22]W $fieldinfo[$fieldoffset+11]d$fieldinfo[$fieldoffset+13].$fieldinfo[$fieldoffset+15]N\' | cs2cs +proj=latlong +datum=NAD83 +to +proj=lcc +datum=NAD27 +units=ft +lon_0=-82.5 +lat_0=39.666666667 +lat_1=40.433333333 +lat_2=41.433333333 +x_0=609601.21920 +y_0=0 > /tmp/coordinates");
      #
      # Open the coordinate file.
      open(COORDS,"</tmp/coordinates");
      #
      # Read the coordinates in.
      my @coords = split('\t', <COORDS>);
      my @northelev = split(' ',$coords[1]);
      #
      # Print out the coordinates to see if we have this right.
      #print "Easting = $coords[0], Northing = $northelev[0], Elevation = $fieldinfo[$fieldoffset+16]\n";
      #
      # Close the coordinate file.
      close COORDS;
      #
      # Set the x & y for the point object.
      $pnt->{x} = $coords[0];
      $pnt->{y} = $northelev[0];
      #
      # Add the point to the shapefile.
      $tcounts->addPoint($pnt);
    }
    #
    # Print the unpack string to see if we got this right.
    #print "UnPack String = $unpackstr\n";
    #
    # Add the next fields standard 10 parts to the unpack string.
    $unpackstr = $unpackstr . " A19 a10 n n a N N N a14 n";
    #
    # Record where this field started at.
    $offsets[$fld] = $fieldoffset;
    #
    # Set the field offset to include the now completed field.
    $fieldoffset = $fieldoffset + $parts[$fieldinfo[$fieldoffset+2]];
  }
  #
  # Does the dbf need created and is this the first record.
  if ( ( $dbfreccnt == -1 ) && ( $recrd == 1 ) ) {
    #

    # Set the record count to 0.
    $dbfreccnt = 0;
    #
    # How many fields are there.
    my $fldcnt = scalar(@offsets);
    #
    # Initialize the field names, type, length, & decimal strings to blank.
    my $fldnames = '';
    my $fldtypes = '';
    my $fldlenth = '';
    my $flddecml = '';
    #
    # Loop through each field and concatenate the name, type, length, & decimal together.
    for ($iname=0; $iname<$fldcnt; $iname++) {
      #
      # Escape out any unprintable characters in the field name.
      $fieldinfo[$offsets[$iname]+0] = uc($fieldinfo[$offsets[$iname]+0]);
      $fieldinfo[$offsets[$iname]+0] =~ s/[^\x41-\x5A]//g;
      #
      # The field type comes in strange sometimes so this should truncate it
      #   so it only contains values of 0-10.
      #   (binary/unpack guru applications now being accepted).
      # Note: Basically this divides by base 16 to move the number 2 decimal
      #       places left then truncates the whole number and multiplies by
      #       base 16 to move the decimal 2 places right.
      $fieldinfo[$offsets[$iname]+2] = 256 * ( ( $fieldinfo[$offsets[$iname]+2] / 256 ) - ( int( $fieldinfo[$offsets[$iname]+2] / 256 ) ) );
      #
      # Okay, the same thing happens with the range.
      $fieldinfo[$offsets[$iname]+3] = 256 * ( ( $fieldinfo[$offsets[$iname]+3] / 256 ) - ( int( $fieldinfo[$offsets[$iname]+3] / 256 ) ) );
      #
      # Concatenate the field name.
      $fldnames = $fldnames . ' "' . $fieldinfo[$offsets[$iname]+0] . '"';
      #
      # Concatenate the field types.
      $fldtypes = $fldtypes . ' "' . $dbfftype[$fieldinfo[$offsets[$iname]+2]] . '"';
      #
      # Concatenate the field lengths.
      $fldlenth = $fldlenth . ' "' . $dbffsize[$fieldinfo[$offsets[$iname]+2]] . '"';
      #
      # Concatenate the field decimals.
      #   All undef right now.
      $flddecml = $flddecml . ' "undef"';
      #
      # If this is not the last field throw in a comma.
      if ( $iname != ( $fldcnt - 1 ) ) {
        $fldnames = $fldnames . ',';
        $fldtypes = $fldtypes . ',';
        $fldlenth = $fldlenth . ',';
        $flddecml = $flddecml . ',';
      }
    }
    #
    # Add the fields for the record number and error flag.
    $fldnames = $fldnames . ', "RECORD", "ERRFLAG"';
    $fldtypes = $fldtypes . ', "N", "N"';
    $fldlenth = $fldlenth . ', "6", "2"';
    $flddecml = $flddecml . ', "undef", "undef"';
    #
    # Create the xbase call.
    my $xbcall = 'XBase->create(name => "' . $sfile . '.dbf", field_names => [' . $fldnames . ' ], field_types => [' . $fldtypes . ' ], field_lengths => [' . $fldlenth . ' ], field_decimals => [' . $flddecml . ' ]) or die XBase->errstr;';
    #
    # Print out the create line to see if we got this right.
    #print "Field Names = $fldnames\nField Types = $fldtypes\nField Sizes = $fldlenth\nField Decimals = $flddecml\n";
    #
    # Create the dbf file.
    $dbh = eval($xbcall);
  }
  #
  # Add the data for this pdb record to the dbf file.
  #
  # Start the xbase add-record call.
  my $xbadd = '$dbh->set_record($dbfreccnt,';
  #
  # How many fields are there.
  my $fldcnt = scalar(@offsets);
  #
  # Go through each field and concatenate the values together.
  for ($iname=0; $iname<$fldcnt; $iname++) {
    #
    # The field type comes in strange sometimes so this should truncate it
    #   so it only contains values of 0-10.
    #   (binary/unpack guru applications now being accepted).
    # Note: Basically this divides by base 16 to move the number 2 decimal
    #       places left then truncates the whole number and multiplies by

    #       base 16 to move the decimal 2 places right.
    $fieldinfo[$offsets[$iname]+2] = 256 * ( ( $fieldinfo[$offsets[$iname]+2] / 256 ) - ( int( $fieldinfo[$offsets[$iname]+2] / 256 ) ) );
    #
    # Is this a number type record.
    if ( $fieldinfo[$offsets[$iname]+2] == 7 ) {
      $xbadd = $xbadd . $fieldinfo[$offsets[$iname]+10];
    }
    #
    # Is this a location type record.
    if ( $fieldinfo[$offsets[$iname]+2] == 4 ) {
      $xbadd = $xbadd . '"' . $fieldinfo[$offsets[$iname]+18] . 'd' . $fieldinfo[$offsets[$iname]+20] . '.' . $fieldinfo[$offsets[$iname]+22] . 'W,' . $fieldinfo[$offsets[$iname]+11] . 'd' . $fieldinfo[$offsets[$iname]+13] . '.' . $fieldinfo[$offsets[$iname]+15] . 'N,' . $fieldinfo[$offsets[$iname]+16] . '"';
    }
    #
        # If this is not the last field throw in a comma.
    if ( $iname != ( $fldcnt - 1 ) ) {
      $xbadd = $xbadd . ',';

    }
  }
  #
  # Add the closer to the end of the xbase add-record call.
  $xbadd = $xbadd . ', ' . $dbfreccnt . ', 0);';
  #
  # Print the xbase add-record line to see if we got this right.
  #print "$xbadd\n";
  #
  # Add the record to the dbf file.
  eval($xbadd);
  #
  # Increment the record counter.
  $dbfreccnt = $dbfreccnt + 1;
}
#
# Close the new shapefile.
undef $tcounts;
#
# Close the dbf handle/file.
undef $dbh;
#
# Get rid of temporary shapefiles if needed.
if ( -e "thistemp.shp" ) {
  unlink "thistemp.shp";
  unlink "thistemp.shx";
  unlink "thistemp.dbf";
}

back to PerlMapScript

Last modified 13 years ago Last modified on Jan 29, 2009, 6:49:43 AM