[GRASS-dev] porting r.in.wms to python

[Glynn Clements]

Hamish wrote:

> semi-related-- does anyone know of a libgis function that will return
> the current projection+coordinate's convergence angle WRT true north?
> A method to find that with the `proj` program is given here:
>  http://grass.osgeo.org/wiki/Ps.map_scripts#Creating_a_fancy_North_Arrow

The -V option is handled by vprocess() in proj.c (in PROJ); after
converting to/from lat/lon, it basically just calls pj_factors() (in
pj_factors.c) then prints out all of the values.

The following code (from r.proj) will get you the "struct pj_info" for
the current location.

    struct pj_info oproj;

    /* Get projection info for output mapset */
    if ((out_proj_info = G_get_projinfo()) == NULL)
	G_fatal_error(_("Unable to get projection info of output raster map"));

    if ((out_unit_info = G_get_projunits()) == NULL)
	G_fatal_error(_("Unable to get projection units of output raster map"));

    if (pj_get_kv(&oproj, out_proj_info, out_unit_info) < 0)
	G_fatal_error(_("Unable to get projection key values of output raster map"));

The "pj" field (i.e. oproj.pj) is the projPJ* which the PROJ library
wants. You should then be able to use pj_factors() directly to get the
convergence angle, e.g.:

	LP lp = { <lat, long> };
	struct FACTORS fact;
	pj_factors(lp, oproj.pj, 0.0, &fac);
	convergence = fac.conv;

-- 
Glynn Clements <glynn at gclements.plus.com>