[pd]:July 8, 2018. This binding to proj.4 uses deprecated functions. [proj] is a replacement binding to the new API using [swig].
[pd]: This a binding to the [http://proj4.org/%|%proj4%|%] projection library. See [https://githubproj4.comrg/OSGdevelo/pment/roj.4eference/wiki/PdeprojAPI ecated.html#example%|% Proj API reference Example %|%]
This implementation requires [http://andreas-kupries.github.io/critcl/%|%critcl%|%]
<<discussion>> 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
======
<<discussion>>
**Compile**
======
critcl -pkg proj4tcl.tcl
======
**Use**
======
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]]]https://www.google.com.au/?gws_rd=ssl
# same, using new instead of create
======
**Using Vertical Datums **
[https://github.com/OSGeo/proj.4/wiki/VerticalDatums%|%Vertical Datums%|%] <<br>>
[http://vdatum.noaa.gov/docs/gtx_info.html%|%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
======
<<discussion>> 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
======
<<discussion>>
**Windows**
Place proj_api.h and libproj-0.dll in the same directory as proj4tcl.tcl.<<br>> 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
======
<<categories>>Geography | Critcl | GIS