[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" [http://www.openmap.org]. The format for shape files is described in an ESRI white paper [http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf] '''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 [list .c scale all %x %y $sqrt2 $sqrt2] bind .c [list .c scale all %x %y $sqrt05 $sqrt05] bind .c [list .c scale all %x %y $sqrt05 $sqrt05] ---- I've refactored [joheid]s 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: * There is also the tclshp extension at http://geology.usgs.gov/tools/metadata/tclshp, which is a wrapper around the free shapelib at http://shapelib.maptools.org. [joheid] * gpsmanshp at http://www.ncc.up.pt/gpsmanshp/ that "provides the means for creating and reading files in the ESRI Shapefile format for keeping 2 or 3 dimensional points and polylines". It is also based on shapelib. * Mapscript/Tcl with MapServer Workbench at http://msworkbench.sourceforge.net/ * the extension [SpatiaLite] to the database [sqlite] can read and write shapefiles from spatialite databases and can e.g. read a shapefile into a virtual table for convenient SQL querying. (It can, by the way, do the same for csv and txtr tables). ---- See also: * [TclWorld], [GIS], [geodesy] ---- [Category File] - [Category Geography]