Changes between Initial Version and Version 1 of PerlMapScriptExamples35ex8

Timestamp:

Jan 29, 2009, 6:49:07 AM (15 years ago)

Author:

jmckenna

Comment:

--

PerlMapScriptExamples35ex8

@@ +1 @@
+= tcounts.pl =
+The url for the data is http://www.highwayengineer.co.medina.oh.us/StickeDB.pdb .
+{{{
+#!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 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";
+}
+}}}