[GRASS5] Re: What proj library version is best?

[Paul Kelly]

Hello Thierry,
(cc to Kergis, GRASS and PROJ lists and Maciek who was also interested in 
this)

First of all a summary of how I understand projections and co-ordinate systems 
and the C language structures used in the various PROJ derivatives to represent 
these.

'Pure' PROJ performs forward and inverse projections from latitude/longitude 
angles to co-ordinates in the projected system described by the 'PJ' structure. 
A PJ struct (created using the pj_init() function) includes a name (e.g. 
proj=tmerc), ellipsoid parameters (a, es or b etc.) and other 
projection-specific parameters such as latitude of origin, false easting and 
northing etc. Something like (from the original projects.h in GRASS CVS):
         /* base projection data structure */
typedef struct PJconsts {
         XY  (*fwd)(LP, struct PJconsts *);
         LP  (*inv)(XY, struct PJconsts *);
         void (*spc)(LP, struct PJconsts *, struct FACTORS *);
         void (*pfree)(struct PJconsts *);
         const char *descr;
         paralist *params;   /* parameter list */
         int over;   /* over-range flag */
         int geoc;   /* geocentric latitude flag */
         double
                 a,  /* major axis or radius if es==0 */
                 e,  /* eccentricity */
                 es, /* e ^ 2 */
                 ra, /* 1/A */
                 one_es, /* 1 - e^2 */
                 rone_es, /* 1/one_es */
                 lam0, phi0, /* central longitude, latitude */
                 x0, y0, /* easting and northing */
                 k0,     /* general scaling factor */
                 to_meter, fr_meter; /* cartesian scaling */
#ifdef PROJ_PARMS__
PROJ_PARMS__
#endif /* end of optional extensions */
} PJ;

Forward and inverse projection are two specific mathematical operations that 
form part of a co-ordinate system conversion but do not do this in themselves.

On Mon, 26 Apr 2004, Thierry Laronde wrote:
> 
> Hello Paul,
> 
> The original CERL sources included a version of the PROJ.4 library being
> provided supplementary to another derived version used for SCS mapgen
> (mapgen has not been updated to use the more recent (as of 1995) PROJ.4
> and Dave Gerdes decided to provide the last PROJ.4 as an add-on).

This is interesting to me as I had thought the extensions to the PROJ version 
in GRASS (two extra files, get_proj.c and do_proj.c and some additions to 
projects.h) were GRASS-specific but it makes more sense that they were part of 
Mapgen and got integrated into GRASS along with other Soil Conservation Service 
enhancements. It is also consistent with an e-mail I remember reading from 
Gerald Evenden in the PROJ mailing list archive asking was Mapgen still part of 
GRASS... that suggested to me he was interested in Mapgen because it included 
some kind of extension to (some old USGS version of) PROJ.

This extension was to add an API (functions pj_get_kv() and pj_do_proj() 
primarily) for performing proper co-ordinate system conversions rather than 
just forward and inverse projections. A general co-ordinate system conversion 
might involve peforming an inverse projection into latitude/longitude and then 
a forward projection into a different system, at the same time accounting for 
perhaps different units used in the two projected systems. A co-ordinate 
system is described by the pj_info struct, which was added to projects.h:
struct pj_info {
       PJ     *pj;
       double meters;
       int    zone;
       char   proj[100];
};

The string pj_info->proj is used to hold the PROJ name for the projection (e.g. 
tmerc) if it is a projected co-ordinate system, or proj=ll for 
latitude/longitude. This addition is needed because of course 
latitude/longitude can be a co-ordinate system in itself as well as an 
intermediate stage between two projected systems.

pj_info->meters holds a conversion factor for the projected units to meters 
(set by unfact=x in the parameter pairs passed to the initialising function 
pj_get_(string|kv)).

I'm really not sure what the zone is for and haven't seen it used anywhere. 
Finally pj_info->pj is a pointer to a PJ struct (as described above) for the 
projected co-ordinate system. If it is lat/long then this is uninitialised; 
zone and meters are set to 0 and 1.0 respectively. The pj_info struct is 
usually created by pj_get_kv(), pj_get_string() or pj_zero_proj(), but there 
were many examples in the GRASS source code where it was created manually for 
a simple lat/long co-ordinate system.

> The author of PROJ.4 has, in 2003, extract the core libraries and made
> some improvements and fixes in what he calls libproj4.

I am fairly sure if you take an old version of GRASS (before 5.0.0pre4) you 
could simply drop the files from Gerald's libproj4 in alongside the get_proj.c 
and do_proj.c and everything will work fine. I don't think the Mapgen version 
modified anything in PROJ, only added those two files and function prototypes 
and struct definition to projects.h.

The problem with that of course is that it won't do datum transformation, 
which as far as I can make out from discussions on the GRASS mailing lists is 
a requirement for most of the uses people want to put the re-projection 
modules to (i.e. [rsv].proj). There are also other modules in GRASS that use 
PROJ but don't need datum transformation, just conversions to/from latitude 
and longitude angles, e.g. r.sun etc.

> 
> Remotesensing.org has taken (old) PROJ.4 maintainance.

Not just maintenance but the extremely useful addition of general co-ordinate 
system conversions incorporating datum transformation. The interesting thing 
is comparing the approach Frank Warmerdam took in extending PROJ to handle 
general co-ordinate system conversions to that chosen by the authors in the 
SCS. I feel Mapgen kept more to the spirit of the PJ struct by retaining it 
as description of only a projected system, and adding their own 'wrapper' 
struct for incorporating details of the units etc. In remotesensing.org proj 
the PJ struct itself was extended from describing purely a projection, to
describing a complete co-ordinate system. proj=longlat or latlong is used in
place of the Mapgen ll, and to_meter in place of unfact for describing the
unit conversion factor. In fact I think this to_meter may have been in place
in later versions of USGS PROJ than the one GRASS/Mapgen was based on.

The datum transformation parameters are considered an attribute of the
co-ordinate system and if the datum transformation parameters are different
then two co-ordinate systems are considered different, even if the
projection parameters and ellipsoid settings are the same. This can
sometimes result in projecting a point between two identical co-ordinate
systems resulting in an (erroneous) change because they have different
datum parameters. This is now confirmed as a feature of the remotesensing.org
PROJ.4 implementation (see PROJ bug no. 368
http://bugzilla.remotesensing.org/show_bug.cgi?id=368 )

> 
> I tend to want to put libproj4 in place of old PROJ.4, but since this is
> a field you know far better than me I'd like to have your advice.
> libproj4 is the library minus the programs. But I'm not fear anymore of
> writing new things or updating old ones...

Yes well as I said above, I think libproj4 should slot in easily enough to
any version of GRASS prior to 5.0.0pre4, so that certainly applies to the 4.x
version KerGIS is based on. Then in the future perhaps datum transformation
could be implemented with the datum information as attributes of the pj_info
(i.e. co-ordinate system) struct rather than the PJ (i.e. projection) struct.

You (or me or whoever ends up making the change) could possibly even re-use
Frank's datum transformation code from remotesensing.org PROJ, in a
re-arranged format.

> 
> TIA.

Yes, that was well in advance :) Sorry for taking so long to reply but I wanted
to try and make my answer clear and useful to other people and a lot of work
intervened.

> 
> PS: I have read your page about datum conversions for France. If you
> need to read documentations in french and find parts you are
> uncomfortable with, tell me I will translate them for you/answer your
> questions (about the meaning of the french writing of course! Not about
> what the text is talking about ;)

I still have not had any feedback on how useful this file is so I will give
the web address here in case anyone from France is interested:
http://www.stjohnspoint.co.uk/gis/france.htm

Good luck in the continuing efforts with KerGIS. I know there is a lot of
background work in tidying up the code structure and preparing it all for
one day when it will be used by many other developers in diverse GIS projects;
hopefully this will not be too far away.

Paul