Hello list
I needed the Mollweide and Goode projections in GRASS5.0beta2 for some
global maps I am working on, so I added them. It was not too much work,
because the PROJ4 package that everything is based on is already part of
GRASS, but has never been fully implemented in the g.setproj, r.proj,
v.proj or m.proj modules. I added support for the rest of included
projections as well, because it was not that different.
The week-end is over and I have to go back to normal work. But I would
like to release my (incomplete) patch to the interested now, because I do
not have time to work with this on a regular basis and don't know when/if
I will be able to continue. I hope releasing it will help encourange
testing and clarify some points on the work that has to be done to include
full support for all PROJ4 projections in a future GRASS version. There
are still a lot of gotchas in sanity checks and making things to work
flawlesly with other GRASS modules.
The file patch-grass5.0beta2-mh1 does the following:
-enables creating and editing locations in any of the projections of the
proj version that is currently included in grass5.0beta2. With this
the total number of projections supported by GRASS is 76.
-enables reprojection of rasters and vectors with v.proj and r.proj
into the new projections.
WARNING AND DISCLAIMER
This is experimental, incomplete, unofficial and probably buggy code.
I do not take any responsibility for what may happen to your data (or
even to your hardware :v) if you apply this patch. But if you know what
you are doing you will get additional projections to play around with
in GRASS. But be prepared to revert to the official version at its next
release. This patch should only be applied to a clean grass5.0beta2 source
tree. (But changes you have made to files not mentioned in the CHANGES
section does not matter.)
WHERE TO GET IT
You can download the patch from
http://www.ngb.se/~morten/grass/patch-grass5.0beta2-mh1.tar.gz
HOW TO APPLY
1) Untar the patch-grass5.0beta2-mh1.tar.gz file into the top level of
the source directory of grass5.0beta2.
2) Run the command
patch -p1 < patch-grass5.0beta2-mh1
3) Recompile the sources.
'make install' after removing the file in
src/CMD/next_step
(You don't have to clean your source tree. 'Make' will know
which files need to be recompiled.)
CHANGES
The following files have been changed:
src/libes/gis/projection.h
Extended list of supported projections
src/general/g.setproj/local_proto.h
Added prototypes for new functions in get_num.c
src/general/g.setproj/get_num.c
Added new functions to get values for new options
src/general/g.setproj/table.h
Added headers for new projections and options
src/general/g.setproj/table.c
Added new projection and parameter default values
src/general/g.setproj/main.c
Updated code to build PROJ_INFO and added some sanity
checks for some of the new projections.
I have tried to change as little as possible, though some
of the code looks perculiar to me and I would probably have
done things differently. (But what do I know?)
LIMITATIONS AND BUGS
The input map for r.proj _must_ be in a 'unprojected' (LL) or cylindical
projection. The reason is that only in a cylindical projection can
the N,S,E and W values from cellhd be combined and used directly as
corner coorinates (NW NE SE SW), which is what r.proj does. Of course,
trying to project rasters from a spherical projection map in nonsense,
because the corner co-ordinates of the rectangular area may not be part of
the map at all. This does not affect v.proj; v.proj can project vectors using
any projection as input.
The 'north/east must be larger than south/west' restriction prevents
projections and display of partial maps in some areas of conical and
azimuthal projections. It seems to be a very basic assumpion in
GRASS, so I have not even though seriously about by-passing it.
But it means that only one half of a 360-degree polar cap map can be
projected into or displayed from. You _can_ project and display a
complete map (covering the whole default region), but if the
resolution is too high may not have enough memory to do that
projection. r.proj wants to read the whole input map into memory
becore projecting.
Sanity checks and default values of some of the new projections may
be inclomplete or incorrect. I don't know enough about projections to
be able to provide reasonable defaults or checks for all cases.
(Since PROJ modules _do_ error checking and give diagnosis back,
it would be a better solution to just pass input through to a test
exec and have setup fail if PROJ complains.)
The omerc projection needs special handling of input because of the
alternative ways to use parameters. Haven't finished this.
Some projections have fixed parameters (alsk, gs50). g.setproj will not
ask the user for input for these but will write the fixed parameters directly
into PROJ_INFO. PROJ_INFO should always reflect these values and should never
be changed manually. (The PROJ routines don't care but other GRASS
modules may fail).
In global spherical projections v.proj produces vector maps with artifact
lines accoss the map, connecting (i guess) opposite points on the perimeter,
which are actually part of the same line or polygon.
r.proj sometimes produces empty maps with projecting rasters into
some of the spherical projections. Don't know why.
---------------------------------------------------------------------------
I found two BUGS in the original g.setproj:
In src/general/g.setproj/main.c line 180
buf = G_find_key_value("radius", out_proj_keys);
should read
buf = G_find_key_value("a", out_proj_keys);
otherwise, if the ellipsoid is a 'sphere', the radius 'a' will be
reset to 0.0000000000 if g.setproj is run a second time. This bug
is fixed in this patch, but if you do not want to apply the patch
you may at least want to change line 180 manually.
The second bug is that old projection parameters will not be cleared
from PROJ_INFO if the projection is changed with g.setproj, and the
old projection uses parameters that are not used by the new one. I
haven't fixed that. I doubt the usefulness of allowing users to
change the projections of existing locations anyway, because their
data will become unusable.
NOTES
Maximum cartesian N/S and E/W values are of course calculated differently for
different projections. For example, Orthographic global earth has
N/S/E/W_max=earth radius, while Sinusoidal has N/S_max = earth perimeter/4
and E/W_max = eart perimeter. I use my externally installed proj program
to calculate default region values in meters before setting up a new location.
TO DO
It was my intention to update m.proj as well to support the new
projections, but I did not have the time. Code and headers are
very similar to g.setproj, in fact, at least the headers and
table.c could be merged. Maybe I will give it a try later.
But since I have the PROJ4.3.3 package installed externally I have
had no need for m.proj.
Except for trying to get around the limitations mentioned above, there
are some 45 projections in the PROJ4.3.3 package that are not yet
inlcuded in the GRASS sources. But I hesitate to get into this
deeper than necessary for my own purposes.
AUTHOR
Morten Hulden
on linux 2.2.9, libc6, patched RH6 and Mandrake, i686