[GRASS5] More on Projections

[Frank Warmerdam]

GRASS Luminaries,

A little more I need to know about projections:

 o If the proj number in a raster header, or the WIND/DEFAULT_WIND file is
   99 (PROJECTION_OTHER) does that mean that the definition of the projection
   can only be found in the PROJ_INFO file, and that it will be one of the
   types other than state plane, utm, latlong or XY? 

 o I was looking at the r.proj command to see how it gets the projection 
   definition of the source raster.  It appears to depend on the PROJ_INFO 
   file existing (via a call to G_get_projinfo()) even though this file 
   does not seem to exist on most of the sample datasets on the GRASS Sample
   Dataset page.  It seems to me that r.proj should fallback to looking at
   the proj: code of the raster, or the projection definition in the WIND
   file, but as far as I can tell it does not.

   How should a GRASS program safely determine what the projection of a raster
   is?  

For my purposes, it is convenient to have a conventional PROJ.4 string as a 
way of expressing the complete coordinate system of a raster, including 
units conversions.  For this purpose, I have implemented the following
functions:

/***************************************************************************
 *  char *G_get_cell_as_proj4( cell, mapset );
 *
 *  char *cell           Name of cell to query.
 *  char *mapset         Name of mapset to query.
 *
 *  Returns a PROJ.4 definition of the coordinate system.  The returned
 *  string should be freed with G_free() when no longer needed.  XY coordinate
 *  system, or other errors during processing will return in a return value
 *  of NULL. 
 **************************************************************************/

/***************************************************************************
 *  char *G_kv_to_proj4( in_proj_keys, in_units_keys );
 *
 *  struct Key_Value *in_proj_keys
 *                       List of projection keys, as returned by 
 *                       G_get_projinfo().
 *  struct Key_Value *in_units_keys
 *                       List of units keys, as returned by 
 *                       G_get_projunits().
 *
 *  Returns a PROJ.4 definition of the coordinate system.  The returned
 *  string should be free with G_free() when no longer needed.  Note that
 *  this function will return a PROJ.4.4.2 string, which includes support
 *  for the +proj=longlat.
 **************************************************************************/

In retrospect, it seems having G_get_cell_as_proj4() is unnecessary, and that
I should instead have written a G_get_location_as_proj4() since everyone seems
to agree that the projection definition should be identical in all the 
rasters in a location.  Would it be useful to implement something like:

int G_get_location_coordsys( char *location /* can be NULL for current */, 
   	                     Key_Value **proj_info, Key_Value **unit_info );

that would return the projection, and unit definition for a location (in
proj_info, and unit_info) in a form ready to use with pj_get_kv()?  This
function would use the PROJ_INFO and PROJ_UNITS if available, otherwise it
would fallback to the WIND file?  Then programs like r.proj could be updated
to use this function?  

Finally, I am still in great need of examples of locations with "interesting"
PROJ_INFO and PROJ_UNITS files. 

Best regards,

---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerda at home.com
light and sound - activate the windows | http://members.home.com/warmerda
and watch the world go round - Rush    | Geospatial Programmer for Rent

---------------------------------------- 
If you want to unsubscribe from GRASS Development Team mailing list write to:
minordomo at geog.uni-hannover.de with
subject 'unsubscribe grass5'