[GRASS5] Fixing up projection related code...

[Paul Kelly]

Hello everyone
I've been looking at improving and tidying the datum transformation support
for a while now and made a few changes in my local copy of GRASS. I would
just like a little more discussion and / or suggestions and hints before
commiting the changes to CVS.

On Sun, 2 Feb 2003, Frank Warmerdam wrote:

> The most general solution then would likely be maintaining a list of
> datum shift options for any given datum, each having some descriptive
> text that can be offered to the user when setting up their location
> (in g.setproj?).

This is what I have chosen to go for, implementing a new file
('datumtransform.table') which can contain multiple entries for each
datum. There are 4 fields: datum name, transformation parameters in PROJ.4
syntax, "Where Used" and "Comment". The last two fields should be
descriptive and help the user to pick the appropriate parameters.

It doesn't really make sense the way it is partially currently implemented
in GRASS, to look up a 3-parameter transform from the datum.table and
just use it, because as Frank said each datum may have many different
tranformations associated with it and GRASS shouldn't pretend to be able to
choose the best one for the user when really it needs a human making the
decision. E.g. someone might want to use an old or custom set of parameters
simply for compatibility with existing data.

So after selection, the datum transformation parameters will be stored in
the PROJ_INFO file under the new key 'datumparams'.
E.g. 'datumparams: towgs84=a,b,c' will be equivalent to the old
'dx: a' 'dy: b' and 'dz: c'. The pj_get_kv() wrapper function will still
interpret the old style properly but g.setproj will only write new-style
PROJ_INFO files.

Regarding passing ellipsoid and datum names onto PROJ.4, this will only be
done if GRASS doesn't recognise the names from its own ellipse.table and
datum.table files. Otherwise the underlying parameters will be passed on
but not the names. This will probably solve bug 1047 with the problem with
the PROJ_INFO files r.in.gdal automatically creates but I will have a look
at it soon as well.

> A slightly easier approach would be to have a datum name for each
> possible transformation.  This is what they do in MapInfo for instance.
> So you might have nad27-conus, nad27-ntv1, nad27-HPGN, each represening
> a different approximation for shifting to WGS84.  This is likely what
> I would suggest in your case.  The downside is that if a user select nad27
> it may be difficult to automatically determine a list of all the datum
> names that could apply.

Looking back now I realise this method was sort of partially implemented
in GRASS, e.g. there is an 'a-can' datum for Alaska and Canada, separate
from NAD27. These can't really be purged from the table as it would break
compatibility with existing PROJ_INFO files but we should be careful about
adding new datums to the file in the future. There is no problem with
changing entries in the datumtransform.table file though, as this file is
only looked up when creating projection information.

> I do think that import programs that report projection information, or
> offer to setup new locations (like r.in.gdal) should offer the transformation
> used if possible.

I don't really think that programs that import data should try to
reproject it. This is getting too complicated and then there would be too
many modules to try to keep up to date. The conventional GRASS way seems
to be to import the data into a location with the correct projection, then
reproject it to your target location using [rsv].proj.

Currently g.setproj is only interactive. When we create a cmd version it
can have a flag to only change the datum transformation parameters and not
change anything else. There could be a library function to do this as
well, which could maybe be called by the re-projection programs

Some points I would appreciate if people could comment on:
As far as I can make out a datum has a one-to-many relationship with datum
transformations, and one-to-one with ellipsoid and prime meridian. Neither
the datum lists in GRASS nor PROJ.4 include space for an associated prime
meridian so I'm wondering is this correct. I would like to make the GRASS
format as correct as possible. Probably I should investigate how PROJ.4
uses prime meridians but I don't really have that much spare time at the
minute and would appreciate a hint.

If we were to extend the datum.table format we would need to change the
function which reads it.
At present it is like (in src/libes/gis/datum.c)

 if (sscanf(buf, "%99s \"%255[^\"]\" \"%255[^\"]\" \"%255[^\"]\"",
            name, params, where_used, comment) != 4)
 ...
It would seem to be very hard to add extra fields to this without it
causing problems when reading old-style files, but maybe I'm missing
something obvious. Could we just leave out the check for the return value
being 4?

In GRASS 5.1 there can be proper support for co-ordinate systems other
than State Plane:

On Mon, 26 Aug 2002, Eric G. Miller wrote:

> That info should be in the PROJ_INFO file, IMHO.  Anyway, I thought
> we could have a Coordinate_System directory that g.setproj could
> use.  So far, with my state plane files, I have something like:
>
> Coordinate_Systems/
>    Projected/
>       North_American/
>          State_Plane_1927/
>             State_Plane_Alabama_East
>             State_Plane_Alabama_West
>             ...
>             State_Plane_Wyoming_West_Central
>          State_Plane_1983/
>             ...
>
> This kind of thing is easy enough to add to, and g.setproj could
> just use directory traversals until the user selects a file, then
> just copy the file to the new location...

This sounds like a good idea but what are directory traversals and how
should they be implemented in g.setproj? There should be a GUI version as
well to make it look user-friendly and nice as g.setproj is a lot of
people's first encounter with a GRASS command (mine anyway and I spent
days over it) and it is really a terrible piece of software and very
off-putting. Is it all right to use Tcl/Tk in GRASS 5.1 or will there be
a different GUI?

If there is no consensus I can go ahead without the prime meridians; I'm
sure it can probably be added later.

Paul