# 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

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