[GRASS5] Fixing up projection related code...

[Eric G. Miller]

I've been looking at some of the projection related code, and I'm
trying to come up with a course of action to clean it up in places.
Eventually, this should result in:

1.  Having better defined projection files.
2.  Creation of a set of predefined well known coordinate systems
    files.
3.  A rewrite of g.setproj with these ideas in mind:
    -  Types of coordinate system are XY, geographic, or projected
    -  Ask users if they want to use a predefined coordinate
       system (state plane, UTM, others...) or if they want to
       define their own.  Having a collection of predefined
       coordinate systems where the user can just pick it by
       "name", will simplify user's lives as well as insuring the
       info is correctly specified.
    -  Both projected and geographic require datum/ellipsoid info.
       Ask for datum first, ellipsoid parameters can be looked up
       from it.  Otherwise, just get ellipsoid parameters.
    -  Remove state plane from projection list, it'll be covered
       by the predefined systems above (though still need units...)
    -  Allow user defined datums??
    -  Other...?
4.  Have transparent handling for a variety of datum transformations
    incorporated into the library wrappers for proj4.  Most of the
    necessary work for the transformations has been done by Andreas
    Lange, it's a matter of working in the strategies.  Perhaps the
    programs can pass on a preferred stategy when specified by the
    user (for instance, NADCON vs. Bursa Wolf).  Get rid of GRASS
    copy of proj4 and make a dependency on proj4 >= whatever version.
    It should be a fatal error if a datum transformation is impossible.
5.  Miscellaneous:
    -  Specify vertical datums? MSL, ellipsoid height, others... Mostly
       for documentation, since very few can be translated.
    -  Add missing usfoot/usfeet <-> metric conversions. usinch? usmile?
    -  Add prime meridian as datum parameter.
    -  Other?

So, I hacked up a little perl script to convert the state27 and state83
files to PROJ_INFO files with all of the datum/ellipsoid parameters
present (current g.setproj fails on this).  In thinking about this, I
wondered why PROJ_UNITS exists at all.  Silly for a three line file.
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...

Example generated State Plane file:

# 5010:State_Plane_1927_Alaska_Zone_10
name: State_Plane_1927_Alaska_Zone_10
ellps: clark66
y_0: 0
es: 0.00676865799729109
a: 6378206.4
dx: -22
dy: 157
dz: 176
lat_0: 51
f: 294.9786982
lat_1: 53.8333333333333
lat_2: 51.8333333333333
proj: lcc
lon_0: -176
x_0: 914401.828803658

I would add the following for STP 1927:
unit: usfoot
units: usfeet
meters: 0.30480060960121920243

For STP 1983, it's case by case for the above or just meters...

Anyway, some of the fixing up will probably break programs all over
the place.  I don't want to see switch statements like:

  switch (G_projection()) {
     case XY:
        ...
     case UTM:
        ...
     case STP:
        ...
     case LL:
        ...
     default:
        do_other_thing();
  }

I'd rather see something like:

  switch (G_projection_type()) {
     case PROJ_XY:
        ...
     case PROJ_GEOGRAPHIC:
        ...
     case PROJ_PROJECTED:
        ...
     default:
        some_error_func();
  }

The problems with the first are STP is poorly defined, what's so
special about UTM? and you still don't have very good info about
the projection if you really need it -- in which case, you have
to get the projection key/vals anyway...  But, if you only need
to know if you're working with linear units vs. angular units,
then the second method gives you that (or does it?).

Likewise, the WIND files' idea of projection should change to
reflect this.  Currently, the proj: <num> is next to useless
and zone: <num> isn't applicable in most cases.  What are the
WIND files for? Just the region settings, north/south/east/west
and resolutions.  But arguably, it should be possible to know
the XY vs. LL vs. PROJECTED from this file.

Anyway, I wanted to throw this out for discussion.  I'm tired
of reading/responding to people confused about projections in
GRASS (and they are more confusing than need be).

-- 
Eric G. Miller <egm2 at jps.net>