Updated 2017-12-25 22:33:28 by pd

pd: This a binding to the proj4 projection library. See Proj API reference Example

This implementation requires critcl

  proj4tcl.tcl

 # proj4tcl.tcl

 package require critcl 3.1
 package require critcl::class
 critcl::tcl 8.5

 critcl::clibraries -lproj

 critcl::ccode {
     #include <proj_api.h>
 }

 critcl::class::define proj4 {
    
     type projPJ

     constructor {
         if (objc != 1) {
             Tcl_WrongNumArgs (interp, objcskip, objv-objcskip,
                                "initstring");
             goto error;
         }
 
         char *defn = Tcl_GetString(objv[0]);
         instance = pj_init_plus(defn);

         if (instance == NULL) {
             Tcl_AppendResult(interp, "proj definition error: ",
                                 pj_strerrno(pj_errno), NULL);
             goto error;
         }
     } 

     destructor {
         pj_free(instance);
     }

     method to proc {
         char* dest
         Tcl_Obj* point
     } ok {
         int n;
         double x,y,z;
         Tcl_Obj **d,*l;
         Tcl_CmdInfo destinfo;
         projPJ dest_cs;

         if (Tcl_GetCommandInfo(interp, dest, &destinfo) == 0) {
             Tcl_AppendResult(interp,"Command not found: ",
                               dest, NULL);
             return TCL_ERROR;
         }

         dest_cs = (projPJ) destinfo.objClientData;

         if (Tcl_ListObjGetElements(interp, point, &n, &d) != TCL_OK) {
             return TCL_ERROR;
         }
       
         if (n != 3) {
             Tcl_AppendResult(interp,"not 3D point: ", NULL);
             return TCL_ERROR;
         }

         if (Tcl_GetDoubleFromObj(interp, d[0], &x) != TCL_OK
             || Tcl_GetDoubleFromObj(interp, d[1], &y) != TCL_OK 
             || Tcl_GetDoubleFromObj(interp, d[2], &z) != TCL_OK) {
                 return TCL_ERROR;
         }

         if (pj_is_latlong(instance)) {
             x *= DEG_TO_RAD;
             y *= DEG_TO_RAD;
         }

         if (pj_transform(instance,dest_cs,1,1,&x,&y,&z)) {
             Tcl_AppendResult(interp, "proj error: ",
                                 pj_strerrno(pj_errno), NULL);
             return TCL_ERROR;
         }

         if (pj_is_latlong(dest_cs)) {
             x *= RAD_TO_DEG;
             y *= RAD_TO_DEG;
         }

         l = Tcl_NewListObj(0, NULL);
         Tcl_ListObjAppendElement(interp, l, Tcl_NewDoubleObj(x));
         Tcl_ListObjAppendElement(interp, l, Tcl_NewDoubleObj(y));
         Tcl_ListObjAppendElement(interp, l, Tcl_NewDoubleObj(z));
         Tcl_SetObjResult(interp, l);
         return TCL_OK;
     }
 }

 package provide proj4tcl 1

Compile edit

critcl -pkg proj4tcl.tcl

Use edit

lappend auto_path [pwd]/lib
package require proj4tcl

# Create coordinate systems 
proj4 create mga56  "+proj=utm +zone=56 +south +ellps=GRS80 "
proj4 create gda94  "+proj=latlon +axis=neu +ellps=GRS80"

# A 3d point is just a 3-item list
# order is defined by the proj.4 +axis option
set pos [list 502810 6964520 0]

# Convert the point
set result [mga56 to gda94 $pos]

puts [format "%.5f %.5f %.2f" {*}$result]
# lats and lons are in decimal degrees.
# output -> -27.44279 153.02843 0.00
#
# example from proj api page
proj4 create merc "+proj=merc +ellps=clrk66 +lat_ts=33"
proj4 create latlong "+proj=latlong +ellps=clrk66"
puts [format "%.2f %.2f" {*}[latlong to merc [list -16.0 20.25 0.0]]]
# output -> -1495284.21 1920596.79
# 
set merc [proj4 new "+proj=merc +ellps=clrk66 +lat_ts=33"]
set latlong [proj4 new "+proj=latlong +ellps=clrk66"]
puts [format "%.2f %.2f" {*}[$latlong to $merc [list -16.0 20.25 0.0]]]
# same, using new instead of create

Using Vertical Datums  edit

Vertical Datums
Datum Transformation Grid (.GTX) Overview
set env(PROJ_LIB) [pwd]
# need to make/get .gtx file
proj4 create mga56ag09 "+proj=utm +zone=56 +south +ellps=GRS80  +geoidgrids=ausgeoid09.gtx"
set result [mga56ag09 to gda94 $pos]
puts [format "%.5f %.5f %.3f" {*}$result]
# output -> -27.44279 153.02843 41.901

  example create gtx script

proc main {} {
    set infile "AUSGeoid09_GDA94_V1.01_DOV.txt"
    set outfile "ausgeoid09.gtx"
    set lat [expr {-(45.0+59.0/60.0)}]
    set lon 108.0
    set deltalat  [expr {1.0/60.0}]
    set deltalon  [expr {1.0/60.0}]
    set rows  2280
    set cols  3120

    set in [open $infile r]
    puts [gets $in]
    set out [open $outfile wb]
    puts -nonewline $out [binary format QQQQII \
       $lat $lon $deltalat $deltalon $rows $cols]

    for {set i 0} {$i<$rows} {incr  i} {
        seek $out [expr {($rows-$i-1)*$cols*4+40}]
        for {set j 0} {$j<$cols} {incr j} {
            scan [gets $in]  "GEO %f" N
            puts -nonewline $out [binary format R $N]
        }
    }

    close $in 
    close $out
}

main

Windows edit

Place proj_api.h and libproj-0.dll in the same directory as proj4tcl.tcl.
libproj-0.dll becomes part of the package using the preload command.
package require critcl 3.1
package require critcl::class
critcl::tcl 8.5

set dir [pwd]
critcl::cheaders -I$dir
critcl::include proj_api.h
critcl::clibraries -L$dir
critcl::clibraries $dir/libproj-0.dll
critcl::preload $dir/libproj-0


Compiling with mingw32
critcl -target mingw32 -pkg proj4tcl