** proj **
----
PROJ is a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system (CRS) to another. This includes cartographic projections as well as geodetic transformations. PROJ is released under the X/MIT open source license
[https://proj.org/index.html]
[https://proj.org/development/reference/functions.html]
There's now a github page with this code and some docs and links [https://github.com/proj4tcl/]
Includes TEA compatible build for unixlike os's.
proj extension library usage
======
load ./proj.so
# define proj string
set S "+proj=pipeline "
append S "+step +proj=axisswap +order=2,1 "
append S "+step +proj=unitconvert +xy_in=deg +xy_out=rad "
append S "+step +proj=merc +lat_ts=33 +ellps=clrk66"
puts $S
# create proj object from proj string
set P [proj create $S]
puts $P
# lat 20.25 lon -16
set a [list 20.25 -16]
puts $a
# do forward transformation
# lat lon (degrees) to easting northing
set a [proj fwd $P $a ]
puts $a
# do inverse transformation
# easting northing to lat lon (degrees)
set a [proj inv $P $a]
puts $a
proj destroy $P
======
compile
======
cc -shared -fpic -g -Wall -O2 proj.c -o proj.so -ltclstub8.6 -lproj
======
proj.c
======c
#include <tcl.h>
#include <proj.h>
typedef struct ProjState {
Tcl_HashTable hash; /* List projections by name */
int uid; /* Used to generate names */
} ProjState;
int
ProjCmd(ClientData data, Tcl_Interp * interp, int objc,
Tcl_Obj * CONST objv[])
{
ProjState *statePtr = (ProjState *) data;
int new, index;
char name[20];
char *subCmds[] = {
"create", "destroy", "fwd", "inv",
"crs2crs", "norm", NULL
};
enum Proj {
Create, Destroy, Fwd, Inv,
Crs2Crs, Norm
};
if (objc < 2) {
Tcl_WrongNumArgs(interp, 1, objv, "option ?arg ...?");
return TCL_ERROR;
}
if (Tcl_GetIndexFromObj(interp, objv[1], subCmds,
"option", 0, &index) != TCL_OK) {
return TCL_ERROR;
}
switch (index) {
case Create:
{
if (objc != 3) {
Tcl_WrongNumArgs(interp, 1, objv, "option ?arg ...?");
return TCL_ERROR;
}
PJ *P = proj_create(PJ_DEFAULT_CTX, Tcl_GetString(objv[2]));
if (P == NULL) {
Tcl_AppendResult(interp, "proj definition error: ",
proj_errno_string(proj_errno(0)), NULL);
return TCL_ERROR;
}
statePtr->uid++;
sprintf(name, "proj%d", statePtr->uid);
Tcl_HashEntry *entryPtr =
Tcl_CreateHashEntry(&statePtr->hash, name, &new);
Tcl_SetHashValue(entryPtr, (ClientData) P);
Tcl_SetResult(interp, name, TCL_VOLATILE);
return TCL_OK;
}
case Destroy:
{
if (objc != 3) {
Tcl_WrongNumArgs(interp, 1, objv, "option ?arg ...?");
return TCL_ERROR;
}
Tcl_HashEntry *entryPtr =
Tcl_FindHashEntry(&statePtr->hash, Tcl_GetString(objv[2]));
if (entryPtr == NULL) {
Tcl_AppendResult(interp, "Unknown proj: ",
Tcl_GetString(objv[2]), NULL);
return TCL_ERROR;
}
PJ *P = Tcl_GetHashValue(entryPtr);
Tcl_DeleteHashEntry(entryPtr);
proj_destroy(P);
return TCL_OK;
}
case Fwd:
case Inv:
{
if (objc != 4) {
Tcl_WrongNumArgs(interp, 1, objv, "option ?arg ...?");
return TCL_ERROR;
}
Tcl_HashEntry *entryPtr =
Tcl_FindHashEntry(&statePtr->hash, Tcl_GetString(objv[2]));
if (entryPtr == NULL) {
Tcl_AppendResult(interp, "Unknown proj: ",
Tcl_GetString(objv[2]), NULL);
return TCL_ERROR;
}
PJ *P = Tcl_GetHashValue(entryPtr);
Tcl_Obj **d; int n;
if (Tcl_ListObjGetElements(interp, objv[3], &n, &d) != TCL_OK) {
return TCL_ERROR;
}
if (n < 2) {
Tcl_AppendResult(interp,"not a point: ", NULL);
return TCL_ERROR;
}
PJ_COORD a = proj_coord(0.0, 0.0, 0.0, 0.0);
for (int i=0; i<n; i++) {
if (Tcl_GetDoubleFromObj(interp, d[i], &(a.v[i])) != TCL_OK) {
return TCL_ERROR;
}
}
int dir;
if (index == Fwd) {
dir = PJ_FWD;
} else {
dir = PJ_INV;
}
// preferable to make unit conversions explicit in proj string
// if (proj_angular_input(P, dir)) {
// a.lp.lam = proj_torad(a.lp.lam);
// a.lp.phi = proj_torad(a.lp.phi);
// }
a = proj_trans(P, dir, a);
// preferable to make unit conversions explicit in proj string
// if (proj_angular_output(P, dir)) {
// a.lp.lam = proj_todeg(a.lp.lam);
// a.lp.phi = proj_todeg(a.lp.phi);
// }
Tcl_Obj *listPtr = Tcl_NewListObj(0, NULL);
for (int i=0; i<n; i++) {
Tcl_ListObjAppendElement(interp, listPtr,
Tcl_NewDoubleObj(a.v[i]));
}
Tcl_SetObjResult(interp, listPtr);
return TCL_OK;
}
case Crs2Crs:
{
if (objc != 4) {
Tcl_WrongNumArgs(interp, 1, objv, "option ?arg ...?");
return TCL_ERROR;
}
PJ *P = proj_create_crs_to_crs(PJ_DEFAULT_CTX,
Tcl_GetString(objv[2]),
Tcl_GetString(objv[3]), 0);
if (P == NULL) {
Tcl_AppendResult(interp, "proj definition error: ",
proj_errno_string(proj_errno(0)), NULL);
return TCL_ERROR;
}
statePtr->uid++;
sprintf(name, "proj%d", statePtr->uid);
Tcl_HashEntry *entryPtr =
Tcl_CreateHashEntry(&statePtr->hash, name, &new);
Tcl_SetHashValue(entryPtr, (ClientData) P);
Tcl_SetResult(interp, name, TCL_VOLATILE);
return TCL_OK;
}
case Norm:
{
if (objc != 3) {
Tcl_WrongNumArgs(interp, 1, objv, "option ?arg ...?");
return TCL_ERROR;
}
Tcl_HashEntry *entryPtr =
Tcl_FindHashEntry(&statePtr->hash, Tcl_GetString(objv[2]));
if (entryPtr == NULL) {
Tcl_AppendResult(interp, "Unknown proj: ",
Tcl_GetString(objv[2]), NULL);
return TCL_ERROR;
}
PJ *P = Tcl_GetHashValue(entryPtr);
PJ *Q = proj_normalize_for_visualization(PJ_DEFAULT_CTX, P);
statePtr->uid++;
sprintf(name, "proj%d", statePtr->uid);
entryPtr = Tcl_CreateHashEntry(&statePtr->hash, name, &new);
Tcl_SetHashValue(entryPtr, (ClientData) Q);
Tcl_SetResult(interp, name, TCL_VOLATILE);
return TCL_OK;
}
}
return TCL_OK;
}
void
ProjCleanup(ClientData data)
{
ProjState *statePtr = (ProjState *) data;
PJ *P;
Tcl_HashSearch search;
Tcl_HashEntry *entryPtr = Tcl_FirstHashEntry(&statePtr->hash, &search);
while (entryPtr != NULL) {
P = Tcl_GetHashValue(entryPtr);
Tcl_DeleteHashEntry(entryPtr);
proj_destroy(P);
entryPtr = Tcl_FirstHashEntry(&statePtr->hash, &search);
}
ckfree((char *) statePtr);
}
int
Proj_Init(Tcl_Interp * interp)
{
if (Tcl_InitStubs(interp, "8.1", 0) == NULL) {
return TCL_ERROR;
}
if (Tcl_PkgRequire(interp, "Tcl", "8.1", 0) == NULL) {
return TCL_ERROR;
}
ProjState *statePtr = (ProjState *) ckalloc(sizeof(ProjState));
Tcl_InitHashTable(&statePtr->hash, TCL_STRING_KEYS);
statePtr->uid = 0;
Tcl_CreateObjCommand(interp, "proj", ProjCmd, (ClientData) statePtr,
ProjCleanup); Tcl_PkgProvide(interp, "proj", "0.1.0");
return TCL_OK;
}
// vim: cin:sw=4:tw=75
======
----
[swig] based binding for tcl for the proj library using the new API described [https://proj.org/development/reference/index.html%|%here%|%]. A simple binding to the old API using critcl has been moved to [proj4tcl - deprecated].
*** usage***
Here is the c example from [https://proj4.org/development/quickstart.html%|%proj quick start%|%]
======c
#include <stdio.h>
#include <proj.h>
int main (void) {
PJ_CONTEXT *C;
PJ *P;
PJ_COORD a, b;
/* or you may set C=PJ_DEFAULT_CTX if you are sure you will */
/* use PJ objects from only one thread */
C = proj_context_create();
P = proj_create (C, "+proj=utm +zone=32 +ellps=GRS80");
if (0==P)
return puts ("Oops"), 0;
/* a coordinate union representing Copenhagen: 55d N, 12d E */
/* note: PROJ.4 works in radians, hence the proj_torad() calls */
a = proj_coord (proj_torad(12), proj_torad(55), 0, 0);
/* transform to UTM zone 32, then back to geographical */
b = proj_trans (P, PJ_FWD, a);
printf ("easting: %.3f, northing: %.3f\n", b.enu.e, b.enu.n);
b = proj_trans (P, PJ_INV, b);
printf ("longitude: %.9f, latitude: %.9f\n", b.lp.lam, b.lp.phi);
/* Clean up */
proj_destroy (P);
proj_context_destroy (C); /* may be omitted in the single threaded case */
return 0;
}
======
and a translation to tcl
======tcl
lappend auto_path .
package require proj
set C [proj_context_create]
set p "+proj=utm +zone=32 +ellps=GRS80"
set P [proj_create $C $p]
set a [proj_coord [proj_torad 12] [proj_torad 55] 0 0 ]
set b [proj_trans $P $::PJ_FWD $a]
puts [format "easting: %.3f, northing: %.3f" \
[PJ_ENU_e_get [PJ_COORD_enu_get $b]] \
[PJ_ENU_n_get [PJ_COORD_enu_get $b]]]
set c [proj_trans $P $::PJ_INV $b]
puts [format "longitude: %.9f, latitude: %.9f" \
[[$c cget -lp] cget -lam] \
[[$c cget -lp] cget -phi]]
proc proj_list {a} {
set t [$a cget -xyzt]
return [list \
[$t cget -x] \
[$t cget -y] \
[$t cget -z] \
[$t cget -t] ]
}
puts [format "easting: %.3f, northing: %.3f" \
{*}[proj_list $b]]
puts [format "longitude: %.9f, latitude: %.9f" \
{*}[proj_list $c]]
proj_destroy $P
proj_context_destroy $C
======
*** build ***
**** swig wrapper ****
file: proj.i
======
%module proj
%{
#include <proj.h>
%}
%include <proj.h>
======
****linux****
======
swig -tcl -pkgversion 1.0 -I/usr/include proj.i
gcc -shared -fpic -DUSE_TCL_STUBS proj_wrap.c -o proj.so -ltclstub8.6 -lproj
======
**** windows****
using [http://www.msys2.org/%|%msys2/mingw64%|%] static (tested with [http://www.bawt.tcl3d.org/download.html#tclbi%|%bawt%|%])
======
pacboy -S gcc swig proj tcl
swig -tcl -pkgversion 1.0 -I/mingw64/include proj.i
gcc -shared -fpic -DUSE_TCL_STUBS proj_wrap.c -o proj.so -L/mingw64/lib -ltclstub86 -l:libproj.a
======
****both****
file: pkgIndex.tcl
======
package ifneeded proj 1.0 [list load [file join $dir proj.so]]
======
***bugs***
A number of functions are not supported using this simple wrapper namely:
======c
PJ *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv);
int proj_trans_array (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord);
size_t proj_trans_generic (
PJ *P,
PJ_DIRECTION direction,
double *x, size_t sx, size_t nx,
double *y, size_t sy, size_t ny,
double *z, size_t sz, size_t nz,
double *t, size_t st, size_t nt
);
double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coord);
void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION logf);
const PJ_OPERATIONS *proj_list_operations(void);
const PJ_ELLPS *proj_list_ellps(void);
const PJ_UNITS *proj_list_units(void);
const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void);
double proj_dmstor(const char *is, char **rs);
======
<<categories>>GIS