Reading GIS shape files

Arjen Markus (29 january 2004) In an attempt to enhance TclWorld with the ability to read so-called shape files I wrote two scripts, one to deal with binary files and one to read shape files.

The first script is quite general, the second is somewhat limited - it deals with polygons only. And neither gracefully deals with errors ;).

The file that I used is a shape file from "OpenMap" [L1 ].

The format for shape files is described in an ESRI white paper [L2 ]

4Oct2005: Added a "namespace import" command to the second piece; the code doesn't work without it. --Scott Hill


 # binaryfile.tcl --
 #    Library of procedures to easier deal with binary files
 #    The procedures can open the file in binary mode, read
 #    the data in "clusters" and convert them to ordinary lists of
 #    data according to the given format.
 #

 # BinaryFile --
 #    Namespace for the variables and procedures
 #
 namespace eval ::BinaryFile {
    namespace export openBinary readData skipBytes
 }

 # openBinary --
 #    Open the file in binary mode
 # Arguments:
 #    filename   Name of the file to read
 #    mode       Read (r) or write (w) mode (optional, default: read)
 # Result:
 #    Handle to the opened file
 #
 proc ::BinaryFile::openBinary { filename {mode r} } {
    set outfile [open $filename $mode]
    fconfigure $outfile -translation binary
    return $outfile
 }

 # skipBytes --
 #    Skip a number of bytes
 # Arguments:
 #    infile        Handle to the file to be read
 #    skip          Number of bytes to skip
 # Result:
 #    None
 # Side effects:
 #    Skip a given number of bytes in the file
 #
 proc ::BinaryFile::skipBytes { infile skip } {
    read $infile $skip
    #
    # Seek does not seem to work in quite the same way
    # - end-of-file detection
    # seek $infile $skip current
    return
 }

 # readData --
 #    Read the data in the file and transform them according to the
 #    given format
 # Arguments:
 #    infile        Handle to the file to be read
 #    format        Format to use for reading
 #    number        Number of times to apply the format (each will
 #                  correspond to an element in the returned list)
 # Result:
 #    List of data
 # Side effects:
 #    Advance in the file
 #
 proc ::BinaryFile::readData { infile format number } {
    #
    # Determine how many bytes to read
    # Check for errors!
    #
    if { [string length $format] > 1 } {
       #
       # Hm, only for "a"?
       set length [string range $format 1 end]
    } else {
       set pos    [lsearch -exact {i I f d r R q Q} $format]
       set length [lindex {4 4 4 8 4 4 8 8} $pos]
    }

    set firstf [string index $format 0]

    set result {}
    for {set i 0} {$i < $number} {incr i} {
       set input [read $infile $length]
       if { [lsearch {a i I f d} $firstf] >= 0 } {
          binary scan $input $format value
       } else {
          set value [FloatScan $input $firstf]
       }
       lappend result $value
    }
    return $result
 }

 # FloatScan --
 #    Scan a binary string for floating-point numbers not in
 #    native format
 # Arguments:
 #    input         Input string
 #    format        Format to use for reading (r, R, q, Q)
 # Result:
 #    The transformed value
 # Note:
 #    This will be superfluous if the binary command supports
 #    the new formats
 #
 proc ::BinaryFile::FloatScan { input format } {
    #
    # Reverse the bytes?
    #
    set reverse 0
    if { $::tcl_platform(byteOrder) == "bigEndian" } {
       if { $format == "r" || $format == "q" } {
          set reverse 1
       }
    } else {
       if { $format == "R" || $format == "Q" } {
          set reverse 1
       }
    }

    if { $reverse } {
       set bytes [split $input ""]
       set last_byte [expr {[llength $bytes]-1}]
       set reversed ""
       for {set i $last_byte } {$i >= 0 } {incr i -1} {
          append reversed [lindex $bytes $i]
       }
    } else {
       set reversed $input
    }

    if { $format == "r" || $format == "R" } {
       binary scan $reversed "f" value
    } else {
       binary scan $reversed "d" value
    }
    return $value
 }

(I left out the test code: this dealt with a proprietary/legacy type of file)


 # showmap.tcl --
 #    Read an ArcGIS shape file and show the resulting map
 #

 package require Tk

 source "binaryfile.tcl"
 namespace import ::BinaryFile::*
 canvas .c -width 720 -height 360 -bg white
 pack   .c -fill both

 #
 # Open an ArcGIS shape file and read/draw the shapes
 #
 set infile [openBinary "vmap_area_thin.shp"]

 #
 # Position at byte 100 and read the elements one by one
 # Note:
 # End-of-file detection seems trickier than expected,
 # see comments to "skipBytes"
 #
 seek $infile 100 start
 while { ![eof $infile] } {
    skipBytes $infile 8
    if { [eof $infile] } { break }
    set type      [readData $infile i 1]
    set bounds    [readData $infile q 4]
    set numparts  [readData $infile i 1]
    set numpoints [readData $infile i 1]
    set parts     [readData $infile i $numparts]
    set xycoords  [readData $infile q [expr {2*$numpoints}]]

    .c create polygon $xycoords -fill green -outline black
 }
 close $infile

 #
 # Transform the coordinates to pixel values
 #
 .c scale all   0   0 2 -2
 .c move  all 360 180

 #
 # Add the equator, meridians and such
 #
 .c create line   0 180 720 180 -width 2 -fill darkblue
 .c create line   0 150 720 150 -width 1 -fill blue
 .c create line   0 120 720 120 -width 1 -fill blue
 .c create line   0  90 720  90 -width 1 -fill blue
 .c create line   0  60 720  60 -width 1 -fill blue
 .c create line   0  30 720  30 -width 1 -fill blue
 .c create line   0 210 720 210 -width 1 -fill blue
 .c create line   0 240 720 240 -width 1 -fill blue
 .c create line   0 270 720 270 -width 1 -fill blue
 .c create line   0 300 720 300 -width 1 -fill blue
 .c create line   0 330 720 330 -width 1 -fill blue

 .c create line 360   0 360 360 -width 2 -fill darkblue
 .c create line 660   0 660 360 -width 1 -fill blue
 .c create line 600   0 600 360 -width 1 -fill blue
 .c create line 540   0 540 360 -width 1 -fill blue
 .c create line 480   0 480 360 -width 1 -fill blue
 .c create line 420   0 420 360 -width 1 -fill blue
 .c create line 300   0 300 360 -width 1 -fill blue
 .c create line 240   0 240 360 -width 1 -fill blue
 .c create line 180   0 180 360 -width 1 -fill blue
 .c create line 120   0 120 360 -width 1 -fill blue
 .c create line  60   0  60 360 -width 1 -fill blue

 #
 # Zoom in and zoom out via mouse clicks
 #
 set sqrt2  [expr {sqrt(2.0)}]
 set sqrt05 [expr {sqrt(0.5)}]
 bind .c <Button-1> [list .c scale all %x %y $sqrt2  $sqrt2]
 bind .c <Button-2> [list .c scale all %x %y $sqrt05 $sqrt05]
 bind .c <Button-3> [list .c scale all %x %y $sqrt05 $sqrt05]

I've refactored joheids post about tclshp into a list of Tcl extensions for use with shape files. (I hope that's OK) JMM.

Packages that read shape files:


See also: * TclWorld, GIS, geodesy