Christian G. Tveen wrote:
I had't really realized that units had to be meters (I'm very much used
to only use lat/lng even in projected maps.)
I have two further questions: From where did you obtain the precise
values for the boundaries ? (I saw my rounded values written by v.proj)
1 $ grass5 /opt/grass-data/global5/glynn
2 GRASS:~ > g.region n=57:54:28N s=54:24:11N w=7:38:57E e=15:33:24E
3 GRASS:~ > r.mapcalc one=1
4 GRASS:~ > exit
5 $ grass5 /opt/grass-data/merctest/glynn
6 GRASS:~ > g.region s=7000000 n=8000000 w=-500000 e=1000000
7 GRASS:~ > r.proj in=one location=global5
8 GRASS:~ > g.region rast=one
9 GRASS:~ > g.region -p
Line 6 took a bit of trial and error.
Note: I'm not suggesting this as a standard approach.. I was also
trying to figure out your projection parameters (I initially entered
lat_ts=57, which put your region somewhere around Spitzbergen), so I
was repeatedly projecting maps back and forth between locations.
A more general question is how to obtain the initial WIND/DEFAULT_WIND
file. When starting out, I only have the LL boundaries. Do I manually
calculate the Mercator boundaries for those files, or is there a clever
way ? (g.region initially requires a DEFAULT_WIND file, so somehow I
need a boot-strap)
Use m.proj2 to project some sample coordinates. For projections which
transform rectangles to rectangles (e.g. lat/lon -> Mercator),
projecting the opposing corners will suffice.
For projections which will distort the grid, you need to use knowledge
of the projection to select appropriate points; think about what the
projection will do to a rectangle.
E.g. when projecting lat/lon to a typical (northern hemisphere)
conical projection, a rectangle will be converted to a ring segment
(i.e. a "rectangle" in polar coordinates). The minimum projected Y
coordinate would be obtained by projecting the intersection of the
central meridian and the southern parallel; the maximum Y would come
from either the NE or NW corner, and the min/max X from the SE and SW
corners.
One mechanism to automate this is to create a sites list which samples
the source region boundary, convert those to the destination
projection with s.proj, then set the destination region with
"g.region sites=...". This will work for any likely projection (it
won't work in some extreme cases which "fold" the map, but those
aren't particularly common).
Note: there isn't a command to change the default region; you just
have to overwrite the DEFAULT_WIND file with a suitable WIND file.
FWIW, the default region doesn't have any significance beyond
"g.region -d".
--
Glynn Clements <glynn.clements@virgin.net>