[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''
[KPV] 2025-02-07: Fixed bug where multi-part PolyLines were being drawn as one part.
----
======
# 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::*
proc SplitParts {xycoords parts} {
# Splits xycoords into multiple lists based off indices in parts
if {[llength $parts] == 1} { return [list $xycoords] }
set result {}
set first NONE
set parts2 [lmap x $parts { expr {2 * $x} }]
foreach idx [concat $parts2 end] {
if {$first ne "NONE"} {
if {$idx eq "end"} {
set segment [lrange $xycoords $first end]
} else {
set segment [lrange $xycoords $first $idx-1]
}
lappend result $segment
}
set first $idx
}
return $result
}
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}]]
foreach xy [SplitParts $xycoords $parts] {
.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 [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).
* https://github.com/anoved/Shapetcl%|%Shapetcl%|%, a work-in-progress C extension that provides general-purpose read/write access to ESRI shapefiles. A draft of the API documentation is available https://github.com/anoved/Shapetcl/wiki/Documentation%|%here%|%.
----
See also:
* [TclWorld], [GIS], [geodesy]
<<categories>> File | Geography