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 [L1 ] [L2 ]

There's now a github page with this code and some docs and links [L3 ] 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

#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", "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 here . A simple binding to the old API using critcl has been moved to proj4tcl - deprecated.

usage

Here is the c example from proj quick start

#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

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 msys2/mingw64 static (tested with 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:

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);

GIS