tclshapefile

Difference between version 10 and 12 - Previous - Next
Tclshapefile is a tcl-only package to read esri shapefiles (well at least points, polylines and polygones). It reads both the geometry (.shp) files and the attribute (.dbf) file. Data are returned as [dict]s.
It uses the specific reader components shpreader and [dbfreader]. 
Alternative: [Shapetcl] is a C extension for Tcl based on the Shapefile C Library.
---- ====== # $Id: tclshapefile.tcl,v 1.1 2012/01/02 16:18:03 joheid Exp $ # # package to read esri shapefiles # uses shpreader and dbfreader as components # (C) Joachim Heidemeier 2011 # published under the same license als Tcl # # methods # readShape # readCoords -> delegated # readAttributes -> delegated # getCoords -> delegated # getAttributes -> delegated # dbfInfo -> delegatd # shpInfo -> delegated # getRecord # readProjection # package require snit # specific reader packages package require shpreader.tcl package require dbfreader.tcl # # # # # # package provide shapefile 0.5 snit::type shapefile { # # components # component dbfReader delegate option -encoding to dbfReader delegate method attributes to dbfReader delegate method readAttributes to dbfReader as readData delegate method getAttributes to dbfReader as data delegate method dbfInfo to dbfReader as dbfInfo # component shpReader delegate option -bbox to shpReader delegate option -shapetype to shpReader delegate option -nrecords to shpReader delegate method readCoords to shpReader as readData delegate method getCoords to shpReader as coords delegate method shpInfo to shpReader # # options # # option -filename -validatemethod CheckFile -configuremethod ShapeFileName # # for use as component # option -partof # option -projection none # # variablen # variable status ;# interneral status set status(shpcheck) 0 set status(dbfcheck) 0 set status(projection) 0 set status(projectionfile) {} # # validatemethods # # checks whether a shp, shx and dbf and proj file is available # error if no shp & dbf # method CheckFile {option value} { set fn [file rootname $value] if {!([file exists $fn.shp] && [file exists $fn.shx] && [file exists $fn.dbf])} { error "$fn ist kein gültiger shapefile" } else { set status(shpcheck) 1 set status(dbfcheck) 1 } if {[file exists $fn.prj]} { set status(projection) 1 } } # # configuremethods # method ShapeFileName {option value} { set f [file rootname $value] set options(-filename) $f set status(projectionfile) ${f}.prj $dbfReader configure -filename ${f}.dbf $shpReader configure -filename ${f}.shp return } # # constructor {args} { install dbfReader using dbasefile %AUTO% -partof $self install shpReader using shpfile %AUTO% -partof $self $self configurelist $args } destructor { catch {$dbRreader destroy} catch {$shpReader destroy} } # # read the proj file # public methods method readProjection {} { if {$status(projection)} { set option(-projection) [read [open $status(projectionfile) r]] } } # # # opens shp, dbf and proj, set options # reads all data # method readShape {} { $shpReader readData $dbfReader readData if {[$shpReader cget -nrecords] != [lindex [$dbfReader dbfInfo] 0]} { error {different number of records in shp and dbf File} } $self readProjection } method getRow {recnr} { set x [dict create] if {[string is integer $recnr] && $recnr <= [$shpReader cget -nrecords]} { dict append x coords [dict get [$shpReader coords $recnr] $recnr] dict append x attributes [dict get [$dbfReader data $recnr] $recnr] } else {error "not existing record requested"} return $x } } # # if {0} { shapefile s1 -filename [tk_getOpenFile] s1 readShape set c [s1 getCoords 1] set a [s1 getAttributes 1] set x [s1 getRow 1] set t [s1 cget -shapetype] } ====== ---- Now the shp-file reader: ====== # $Id: shpreader.tcl,v 1.3 2012/01/02 16:17:40 joheid Exp $ # (C) Joachim Heidemeier # published under the same licence as Tcl # pure tcl snit object to read the files with the geometries in shapefiles (.shp files) # public methods # readData -- read and decode the file # coords -- outputs the coordinates of one or all geometric objects # shpInfo -- give information on internal details (mostly for inspection purposes) # shapetype, bounding box and number of geometric objects are available as (read-only) # obtions. # if {0} { The information on shapefile structure is taken from the ESRI shapefile documentation (http://www.esri.com/library/whitepapers/pdfs/shapefile.pdf). In the tables the first columns including ByteOrder/Content are taken from the documentation. A a shapefile consists of a 100 Byte long Header and subsequent records made up from an 8 byte record header and the record content. The record content structure is shapetype dependent. 1. Definition of Shape File Header (100 Bytes long), {II5Iiiqqqqqqqq} Position Field Value Type ByteOrder FormatSpecifier TclVariable ===================================================================================================== Byte 0 FileCode 9994 Integer Big I info(fc) Byte 4 Unused 0 Integer Big { discarded Byte 8 Unused 0 Integer Big discarded Byte 12 Unused 0 Integer Big I5 discarded Byte 16 Unused 0 Integer Big discarded Byte 20 Unused 0 Integer Big } discarded Byte 24 FileLength FileLength Integer Big I info(FLength),in Bytes Byte 28 Version 1000 Integer Little i discarded Byte 32 ShapeType ShapeType Integer Little i option(-shapetype) Byte 36 BoundingBox Xmin Double Little q { Byte 44 BoundingBox Ymin Double Little q Byte 52 BoundingBox Xmax Double Little q option(-bbox) Byte 60 BoundingBox Ymax Double Little q as dict, key coordname Byte 68* BoundingBox Zmin Double Little q { Byte 76* BoundingBox Zmax Double Little q if required, Byte 84* BoundingBox Mmin Double Little q else discarded Byte 92* BoundingBox Mmax Double Little q }} ===================================================================================================== * Unused, with value 0.0, if not Measured or Z type The value for file length is the total length of the file in 16-bit words (including the fifty 16-bit words that make up the header). Available shape types Value ShapeType Content supported =================================================== 0 NullShape Null y 1 Point P-Struct y 3 PolyLine PL-Struct y 5 Polygon PG-Struct y 8 MultiPoint 11 PointZ 13 PolyLineZ 15 PolygonZ 18 MultiPointZ 21 PointM 23 PolyLineM 25 PolygonM 28 MultiPointM 31 MultiPatch ====================================================== 2. Record Header fixed (8 bytes long) Position Field Value Type ByteOrder FormatSpecifier TclVariable ================================================================================================================================ Byte 0 RecordNumber RecordNumber Integer Big I IInfo(RecNR), Key für coords Dict Byte 4 ContentLength ContentLength Integer Big I info(RecNR,Rlength), in Bytes ================================================================================================================================= The content length for a record is the length of the record contents section measured in 16-bit words. Each record, therefore, contributes (4 + content length) 16-bit words toward the total length of the file, as stored at Byte 24 in the file header. 3. Record Content (variable length) ShapeType(4Bytes) + geometric data according to type returned Tcl Variable coords is a dict with RecordCount entries, key is RecNr value is list of X Y coordinates 3.1 Point-Struct Position Field Value Type Number ByteOrder FormatSpecifier TclVariable ============================================================================================================== Byte 0 ShapeType 1 Integer 1 Little i internal check Byte 4 X X Double 1 Little q coords Byte 12 Y Y Double 1 Little q coords ============================================================================================================== 3.2 Polyline-Struct consists of parts Position Field Value Type Number ByteOrder Remark TclVariable =============================================================================================================== Byte 0 ShapeType 3 Integer 1 Little internal check Byte 4 BBox BBox Double 4 Little xmin,ymin,xmax,ymax info(bbox,recNr) Byte 36 NumParts NumParts Integer 1 Little info(numParts,recNr) Byte 40 NumPoints NumPoints Integer 1 Little for all parts info(numPoints,recNr) Byte 44 Parts Parts Integer NumParts Little offset to first point of part, 0_indexed Byte X Points Points Point NumPoints Little no part delimiter Note: X = 44 + 4 * NumParts ================================================================================================================ 3.3 Polygon-Struct consists of rings, structure analogue to PL. Position Field Value Type Number ByteOrder Remark TclVariable ========================================================================================================== Byte 0 ShapeType 5 Integer 1 Little internal check Byte 4 BBox BBox Double 4 Little xmin,ymin,xmax,ymax Byte 36 NumParts NumParts Integer 1 Little Byte 40 NumPoints NumPoints Integer 1 Little for all rings Byte 44 Parts Parts Integer NumParts Little offset to first point of part, 0_indexed Byte X Points Points Point NumPoints Little no part delimiter Note: X = 44 + 4 * NumParts ============================================================================================================= } # # # package require snit # package provide shpreader 0.5 snit::type shpfile { variable TCL_OK 0 variable TCL_CONTINUE 3 # # # daten as list of lists variable coords variable info variable fh variable status ;# interner Status variable shapeTypes ;# gueltige shapes mit Nr / Namenmapping set status(shp_ok) 0 set status(shp_open) 0 set status(eof) 0 # options option -filename -configuremethod ShpFileName # ohne extension, die wird in configuremethod weggeworfen option -bbox -readonly yes ;# bounding Box wird aus shapefile gelesen # normally used as component option -partof # shapetype option -shapetype -default null -readonly yes option -nrecords -readonly yes # # configuremethods # method ShpFileName {option value} { # set options(-filename) "[file rootname $value].shp" return } # constructor {args} { # set status(read) 0 array set shapeTypes [list\ 0 Null\ 1 Point\ 8 MultiPoint\ 3 PolyLine\ 5 Polygon\ 21 PointM\ 28 MultiPointM\ 23 PolyLineM\ 25 PolygonM\ 11 PointZ\ 18 MultiPointZ\ 13 PolyLineZ\ 15 PolygonZ\ 31 MultiPatch] $self configurelist $args } # # reads and evalutes the header # stores relevant data in the info or options variable # method ReadHeader {} { # formatstring for binary scan accoring to table 1 set formatstring {II5Iiiqqqqqqqq} # set rec [read $fh 100] set status(eof) [eof $fh] if {!$status(eof)} { binary scan $rec $formatstring info(fc) tmp info(fLength) vers info(type) xm ym XM YM zm ZM mm MM # if {$info(fc) != 9994} { error {file is not a shp-file} } else { set info(fLength) [expr {2 * $info(fLength)}] set options(-shapetype) $shapeTypes($info(type)) set options(-bbox) [dict create XMin $xm YMin $ym XMax $XM YMax $YM] if {$info(type) > 8} { if {$info(type) < 21} { dict append options(-bbox) ZMin $zm dict append options(-bbox) ZMax $ZM } else { dict append options(-bbox) MMin mm dict append options(-bbox) MMax MM } } return $TCL_OK } } else { return $TCL_CONTINUE } } # # open File or raise error # method OpenShpFile {} { if {[catch {open $options(-filename)} res]} { return -code error $res } set fh $res set status(shp_open) 1 fconfigure $fh -translation binary return $TCL_OK } # # decode a point record # proc getPoint {rcont} { if {[string length $rcont != 20]} {error "wrong contentlength in getPoint [string length $rcont] instead of 20"} binary scan $rcont iq2 t c if {$t != 1} {error "wrong shapetype in getPoint $t instead 1"} return [dict create p $c] } # # Polyline and Polygone need the same algorithm # returns a dict over parts # proc getPoly {selfns type recnr rcont} { set pp 0 set offset 0 set tmp [dict create] # get type and check, we assume one type per file binary scan $rcont i t if {$t != $type} {error "wrong shapetype in getPoly $t instead $type"} set offset 4 ;#start of data section # get record info binary scan $rcont @${offset}q4ii bbox npart npoint set info(bbox,$recnr) $bbox set info(numParts,$recnr) $npart set info(numPoints,$recnr) $npoint set offset 44 ;#start of variable section # get part-offsets set formatstring "@${offset}i$npart" binary scan $rcont $formatstring ofs set offset [expr {44 + 4* $npart }] # set end of last part lappend ofs $npoint # build formatstring and varlist for scan set formatstring "@$offset" set varlist [list] set keylist [list] for {set i 0} {$i < $npart} {incr i} { set start [lindex $ofs $i] set end [expr {[lindex $ofs $i+1] -1}] set np [expr {$end - $start} +1] set nfloat [expr {$np * 2}] append formatstring "q$nfloat" set pp [expr {$pp + $np}] lappend keylist part${i} } set varlist $keylist # # the conversion # binary scan $rcont $formatstring {*}$varlist # # the return value # foreach key $keylist var $varlist { dict append tmp $key [set $var] } if {$pp != $npoint} {error "wrong number of points"} return $tmp } # # method GetNextRecord # reads one record # assign it to dict # method GetNextRecord {} { # first the RecordHeader set headerFormat II set header [read $fh 8] set status(eof) [eof $fh] if {!$status(eof)} { binary scan $header $headerFormat recNr contentLength set contentLengthBytes [expr {$contentLength *2}] } else { return $TCL_CONTINUE } set content [read $fh $contentLengthBytes] set status(eof) [eof $fh] if {!$status(eof)} { switch -exact -- $info(type) { 0 {dict append coords $recNr {}} 1 { dict append coords "$recNr" [getPoint $content] } 3 { dict append coords "$recNr" [getPoly $selfns $info(type) $recNr $content] } 5 { dict append coords "$recNr" [getPoly $selfns $info(type) $recNr $content] } default {error "unsupported shape type in GetNextRecord"} } set info(recnr) $recNr return $TCL_OK } else { return $TCL_CONTINUE } } # # public methods # read the file # method readData {} { if {!$status(read)} { $self OpenShpFile $self ReadHeader set retval [$self GetNextRecord] while {$retval == $TCL_OK} { set retval [$self GetNextRecord] } set options(-nrecords) $info(recnr) set status(read) 1 if {[catch {close $fh} msg]} { unset fh return -code error -errorcode $errorCode $msh } return $TCL_OK } else { error {file already read} } } # # retrieve the coords # of one record or all # argument # # method coords {{recnr all}} { if {"$recnr" eq "all"} { return $coords } elseif {[string is integer $recnr] && $recnr <= $info(recnr)} { return [dict filter $coords key $recnr] } else {error "Record $recnr not available"} } # # return internal infos from info array # what is optional pattern # method shpInfo {{what all}} { if {"$what" eq "all"} { return [array get info] } else { set s [array get info $what] if {[llenth $s]} { return $s } else { return "info $what not available" } } } } ====== ---- '''[arjen] - 2021-09-15 06:08:48''' I tried to use this package to read a shapefile and use the coordinates to visualise some modelling results. Unfortunately, the shapefile in question was not supported, due to unknown item types. However, it so happened that the accompanying dbf file (which this package can also read) contained enough information for me to do the visualisation. ---- '''[arjen] - 2021-09-16 06:55:16''' I had the opportunity to try and read yet another sahapefile: this time all elements it contained are supported, but apparently there is no projection information as I got the message: ====== can't read "status(projection)": no such element in array while executing "if {$status(projection)} { set option(-projection) [read [open $status(projectionfile) r]] }" (procedure "::shapefile::Snit_methodreadProjection" line 6) invoked from within "$self readProjection" (procedure "::shapefile::Snit_methodreadShape" line 11) invoked from within "s1 readShape" (file "tclshape.tcl" line 141) ====== As I do not use any GIS system (more for lack of experience than for a more valid reason), but am confronted from time to time with shapefiles, I wonder if I should try and update this code ... ;) ---- '''Harm Olthof - 2021-09-16'''
Note that for the Dutch governemental context the use of shp-files are more or lesse forbidden since the modern alternative for the ancient shp format, the https://www.geopackage.org/%|%Geopackage%|% is on the "explain-or-comply" list of Forumstandaardisatie. Since a Geopackage is basically an [SQLite] database with spatial extensions, this may be an easier way out: either install the Open Source desktop GIS-application https://www.qgis.org%|%QGIS%|% and drag and drop the .shp in the canvas, (in case of a missing .prj you'll get a popup where you can provide the projection information) and then right-click the layer and export it to a Geopackage or use the command line to do this with the gdal library (for which, a long time ago existed tcl-bindings). See examples https://gdal.org/programs/ogr2ogr.html%|%here%|% Once you genererated the GeoPackage, you can access it with [tdbc] with all it's possibilities.
Note that for the Dutch governemental context the use of shp-files are more or lesse forbidden since the modern alternative for the ancient shp format, the https://www.geopackage.org/%|%Geopackage%|% is on the "comply-or-explain" list of Forumstandaardisatie. Since a Geopackage is basically an [SQLite] database with spatial extensions, this may be an easier way out: either install the Open Source desktop GIS-application https://www.qgis.org%|%QGIS%|% and drag and drop the .shp in the canvas, (in case of a missing .prj you'll get a popup where you can provide the projection information) and then right-click the layer and export it to a Geopackage or use the command line to do this with the gdal library (for which, a long time ago existed tcl-bindings). See examples https://gdal.org/programs/ogr2ogr.html%|%here%|% Once you genererated the GeoPackage, you can access it with [tdbc] with all it's possibilities.
<<categories>> Package | Geography