[GRASS-dev] problems with datum parameters tables

GRASS GRASS developer
1

On Oct 13, 2009, at 3:51 AM, grass-dev-request@lists.osgeo.org wrote:

Date: Tue, 13 Oct 2009 02:18:32 -0700 (PDT)
From: Hamish <hamish_b@yahoo.com>
Subject: Re: [GRASS-dev] problems with datum parameters tables
To: list GRASS developers <grass-dev@lists.osgeo.org>, Michael Barton
       <michael.barton@asu.edu>
Message-ID: <230683.60053.qm@web110011.mail.gq1.yahoo.com>
Content-Type: text/plain; charset=utf-8

Michael Barton wrote:

I'm working on the GUI location
manager and have run into a couple of irregularities in a
couple of the projections parameters files. I'm sending this
to the list instead of keeping it in the location wizard
track ticket because they are a more general issue that will
affect projection/location creation by other means too.

In proj-parms.table, the entry for "GEOS" is incorrect.
There should be 3 subentries for each entry: projection
code, projection description, and required parameters.
However, GEOS is missing the projection description.

GEOS:HEIGH=ask,35800000.0;LON0=ask,20.0

The following entry has the correct format

GINS8:Ginsburg VIII (TsNIIGAiK):LAT0=ask,0.0;LON0=ask,20.0

GEOS is the only incorrect formatted entry in the file.

That's the view-from-a-geostationary-satellite projection. (e.g. Google
Earth app zoomed way out from space) I guess the missing name is a typo.
http://www.remotesensing.org/geotiff/proj_list/geos.html

Need to add this so that we can use it. If it has the wrong number of parameters it won't be parsed correctly.

It's a really neat one to reproject world wide datasets into.
View of the Earth/Mars/Moon from space from any angle and any distance.
FWIW at those distances 100m offset from a datum probably ain't much to
worry about.

Sounds cool.

It also is NOT in the datums.table file.

It shouldn't be. Generally speaking, projections and datums are completely
independent animals and can be mixed at will. Popular convention means that
in practice some datums are always used with some projections, but that's
not to satisfy the math, just history.

I doubt any of the proj-*.table entries will be found in the datum.table
file, except perhaps as a footnote comment.

I meant to say that it is not in the projections file, not the datum.table file. But this brings up a related issue. Except for the missing GEOS projection, it looks to me like all the information in the projections file is also contained in the proj-parms.table file. So it seems like we should be parsing the latter rather than the former for getting projections. What is the point of the projections file anyway if all the same stuff is in the proj-parms.table file?

One difference is that projection codes are in caps in proj-parms.table and in lower case in the projections file. Does this make any difference when making a PROJ4 string?

Also, in the proj-desc.table, several entries are
incorrectly formatted. None of them are used by the
proj-parms.table and so they probably should be removed from
the proj-desc.table file.

IMO fixing is better than removing..

All are missing their description fields that are needed so that they
can be read by a user.

I'm guessing the missing descriptions are due to lack of explanation
of what they do by the person who made the table.

They are:

GUAM:bool:guam:
OALPHA:lat:o_alpha:
OLAT1:lat:o_lat_1:
OLAT2:lat:o_lat_2:
OLON1:lon:o_lon_1:
OLON2:lon:o_lon_2:
OLONC:lon:o_lon_c:

the format looks ok as far as having the correct number of delimiters, the
only issue is that the description is an empty string. Parsing should be
made to deal with that, but descriptions should be found too.

These terms are not unique to GRASS btw, they come from PROJ.4 or prior.
+guam has been around for a long while, you probably have to dig into
the PROJ source code to see what they do exactly (or ask Gerald E).
I've tried to document as many as I know about here:
https://trac.osgeo.org/proj/wiki/GenParms

We will need the human readable description if we ever want to use these parameters. The description is what is presented to the user in order to get him/her to input a value.

Hamish

ps- lest it be forgot, FYI currently we are completely ignoring GRASS's
built-in State Plane support; currently those are only available via the
EPSG files and user knowing+entering custom +proj4 terms.

??? It is in the projections file and hence in the list of projections to choose from in the current location wizard.

Michael

2

Hamish:

> That's the view-from-a-geostationary-satellite projection.
> http://www.remotesensing.org/geotiff/proj_list/geos.html

Michael:

Need to add this so that we can use it. If it has the wrong
number of parameters it won't be parsed correctly.

now fixed in SVN.

I meant to say that it is not in the projections file, not
the datum.table file.

oh, ok. fixed in devbr6 and trunk. no problem in 6.4, it will just
never be looked up.

But this brings up a related issue. Except for the missing GEOS
projection, it looks to me like all the information in the projections
file is also contained in the proj-parms.table file.

(for those playing along at home this is $GISBASE/etc/projections,
see also general/g.setproj/README)

So it seems like we should be parsing the latter rather than the
former for getting projections. What is the point of the projections
file anyway if all the same stuff is in the proj-parms.table
file?

Use the etc/projections file as the master list. the g.setproj
proj-*.table files are lookup tables specifically to aid g.setproj,
but useful for reuse.

One difference is that projection codes are in caps in
proj-parms.table and in lower case in the projections file.
Does this make any difference when making a PROJ4 string?

see the g.setproj README file. they would need to be lower cased
for the string compare/lookup.

>> GUAM:bool:guam:
>> OALPHA:lat:o_alpha:
>> OLAT1:lat:o_lat_1:
>> OLAT2:lat:o_lat_2:
>> OLON1:lon:o_lon_1:
>> OLON2:lon:o_lon_2:
>> OLONC:lon:o_lon_c:

...

We will need the human readable description if we ever want
to use these parameters.
The description is what is presented to the user in order to get
him/her to input a value.

the o_ ones seems to come from the General Oblique Transformation code
in proj4 (PJ_ob_tran.c)

+guam seems to have to do with Azimuthal Equidistant's Guam elliptical
(PJ_aeqd.c)

... http://www.remotesensing.org/geotiff/proj_list/

?

> currently we are completely ignoring GRASS's built-in State Plane
> support;

...

??? It is in the projections file and hence in the list of
projections to choose from in the current location wizard.

State Plane is not a projection itself, it is the name of an index of
projections: there are 1~9 definitions per state. e.g. due to latitude,
geographic extent, and shape, Alaska needs to use a radically different
projection than Michigan, and another totally different projection is
needed for California etc. So 1st the user selects the 1927 or 1983
version of SP defn's, then they select their state, and finally they
choose the County code projection within that.

the county codes are listed in the $GISBASE/etc/FIPS.code file.

There is no +stp proj4 term, it is just an alias within GRASS to tell it
to launch the interactive get_stp.c part of g.setproj.

note in the etc/projection list that ll, utm, and stp are grouped at the
top before any of the "real" projections. (technically LL is unprojected)

I see a bug in g.setproj getting stuck in a loop for the Virgin Islands,
but I suggest to delay work on the SP stuff until after the main datum and
real projection stuff is sorted out. Maybe just have that pop to an error
message for SP saying to use the text version for now.

regards,
Hamish

3

Thanks much for this information. I've got a start on rebuilding the UTM page to show parameters of all projections. I've got the usual formatting hangups (currently can't figure out how to clear the page if you want to look at another projection), but it is moving along.

I agree that we should get this other stuff working and then look into state plane. It may need some kind of enhancement to g.proj.

Michael
____________________
C. Michael Barton
Director, Center for Social Dynamics & Complexity
Professor of Anthropology, School of Human Evolution & Social Change
Arizona State University

Phone: 480-965-6262
Fax: 480-965-7671
www: www.public.asu.edu/~cmbarton, http://csdc.asu.edu

On Oct 13, 2009, at 9:30 PM, Hamish wrote:

Hamish:

That's the view-from-a-geostationary-satellite projection.
http://www.remotesensing.org/geotiff/proj_list/geos.html

Michael:

Need to add this so that we can use it. If it has the wrong
number of parameters it won't be parsed correctly.

now fixed in SVN.

I meant to say that it is not in the projections file, not
the datum.table file.

oh, ok. fixed in devbr6 and trunk. no problem in 6.4, it will just
never be looked up.

But this brings up a related issue. Except for the missing GEOS
projection, it looks to me like all the information in the projections
file is also contained in the proj-parms.table file.

(for those playing along at home this is $GISBASE/etc/projections,
see also general/g.setproj/README)

So it seems like we should be parsing the latter rather than the
former for getting projections. What is the point of the projections
file anyway if all the same stuff is in the proj-parms.table
file?

Use the etc/projections file as the master list. the g.setproj
proj-*.table files are lookup tables specifically to aid g.setproj,
but useful for reuse.

One difference is that projection codes are in caps in
proj-parms.table and in lower case in the projections file.
Does this make any difference when making a PROJ4 string?

see the g.setproj README file. they would need to be lower cased
for the string compare/lookup.

GUAM:bool:guam:
OALPHA:lat:o_alpha:
OLAT1:lat:o_lat_1:
OLAT2:lat:o_lat_2:
OLON1:lon:o_lon_1:
OLON2:lon:o_lon_2:
OLONC:lon:o_lon_c:

...

We will need the human readable description if we ever want
to use these parameters.
The description is what is presented to the user in order to get
him/her to input a value.

the o_ ones seems to come from the General Oblique Transformation code
in proj4 (PJ_ob_tran.c)

+guam seems to have to do with Azimuthal Equidistant's Guam elliptical
(PJ_aeqd.c)

... http://www.remotesensing.org/geotiff/proj_list/

?

currently we are completely ignoring GRASS's built-in State Plane
support;

...

??? It is in the projections file and hence in the list of
projections to choose from in the current location wizard.

State Plane is not a projection itself, it is the name of an index of
projections: there are 1~9 definitions per state. e.g. due to latitude,
geographic extent, and shape, Alaska needs to use a radically different
projection than Michigan, and another totally different projection is
needed for California etc. So 1st the user selects the 1927 or 1983
version of SP defn's, then they select their state, and finally they
choose the County code projection within that.

the county codes are listed in the $GISBASE/etc/FIPS.code file.

There is no +stp proj4 term, it is just an alias within GRASS to tell it
to launch the interactive get_stp.c part of g.setproj.

note in the etc/projection list that ll, utm, and stp are grouped at the
top before any of the "real" projections. (technically LL is unprojected)

I see a bug in g.setproj getting stuck in a loop for the Virgin Islands,
but I suggest to delay work on the SP stuff until after the main datum and
real projection stuff is sorted out. Maybe just have that pop to an error
message for SP saying to use the text version for now.

regards,
Hamish

4

Michael wrote:

I've got a start on rebuilding the UTM page to show
parameters of all projections.

ummmm, could you explain that more? inner UTM terms may be interesting
on the summary page, but beyond that the user just cares about "zone 12
North" + datum choice(s). For UTM the projection will always be TM of
course, and other terms will be fixed based on the zone..

I've got the usual formatting hangups (currently can't figure out
how to clear the page if you want to look at another projection),

ie if the back button is pressed?

but it is moving along.

glad to hear it.

I agree that we should get this other stuff working and
then look into state plane.

there is the option to abandon it and tell users to look it up in the
EPSG file instead, but then they have to know what they are looking for
as County names will not be there.

It may need some kind of enhancement to g.proj.

I think there are enough twists and turns in it that get_stp.c would
have to be rewritten in Python & parse the FIPS.code etc files directly.
Looking at the C code, it doesn't seem so bad.

Hamish

5

On Oct 14, 2009, at 10:22 PM, Hamish wrote:

ummmm, could you explain that more? inner UTM terms may be interesting
on the summary page, but beyond that the user just cares about "zone 12
North" + datum choice(s). For UTM the projection will always be TM of
course, and other terms will be fixed based on the zone..

After selecting the projection in the wizard, you go to a new page. Currently, if the projection is UTM, you get to select zone and hemisphere. If it is not UTM, you don't get to select any parameters. The UTM selection is hard wired into the GUI.

I'm changing this so that the the parameters parsed from proj-parms.table appear for whatever projection is chosen--UTM or otherwise.

Michael

6

Michael wrote:

After selecting the projection in the wizard, you go to a
new page. Currently, if the projection is UTM, you get to
select zone and hemisphere. If it is not UTM, you don't get
to select any parameters. The UTM selection is hard wired
into the GUI.

I'm changing this so that the the parameters parsed from
proj-parms.table appear for whatever projection is
chosen--UTM or otherwise.

LL, UTM, and STP are not real projections so hardwiring them
as exceptions wouldn't be cause for alarm; but ok, if a single
method does the trick then no complaints from me...

cheers,
Hamish

7

Hamish wrote:

> After selecting the projection in the wizard, you go to a
> new page. Currently, if the projection is UTM, you get to
> select zone and hemisphere. If it is not UTM, you don't get
> to select any parameters. The UTM selection is hard wired
> into the GUI.
>
> I'm changing this so that the the parameters parsed from
> proj-parms.table appear for whatever projection is
> chosen--UTM or otherwise.

LL, UTM, and STP are not real projections so hardwiring them
as exceptions wouldn't be cause for alarm; but ok, if a single
method does the trick then no complaints from me...

In particular, I would suspect that using PROJECTION_OTHER for a UTM
projection may confuse either GRASS or the user. E.g. there are many
places where PROJECTION_UTM causes cellhd.zone to become relevant.

OTOH, we might want to consider changing this. There isn't anything
inherently special about UTM; it's not like switching zones is any
different to other projection changes.

Really, there only needs to be 3 PROJECTION_* options: XY, LL and
OTHER. UTM and STP should ultimately be assimilated into OTHER.

--
Glynn Clements <glynn@gclements.plus.com>

8

On Oct 16, 2009, at 10:24 AM, Glynn Clements wrote:

Hamish wrote:

After selecting the projection in the wizard, you go to a
new page. Currently, if the projection is UTM, you get to
select zone and hemisphere. If it is not UTM, you don't get
to select any parameters. The UTM selection is hard wired
into the GUI.

I'm changing this so that the the parameters parsed from
proj-parms.table appear for whatever projection is
chosen--UTM or otherwise.

LL, UTM, and STP are not real projections so hardwiring them
as exceptions wouldn't be cause for alarm; but ok, if a single
method does the trick then no complaints from me...

In particular, I would suspect that using PROJECTION_OTHER for a UTM
projection may confuse either GRASS or the user. E.g. there are many
places where PROJECTION_UTM causes cellhd.zone to become relevant.

OTOH, we might want to consider changing this. There isn't anything
inherently special about UTM; it's not like switching zones is any
different to other projection changes.

Really, there only needs to be 3 PROJECTION_* options: XY, LL and
OTHER. UTM and STP should ultimately be assimilated into OTHER.

-- Glynn Clements <glynn@gclements.plus.com>

I don't know how this works in the underlying C code, this all can be handled in a single method that parses proj-parms.table except for State Plane. I can see why you call latlon a pseudo projection, but not UTM. AFAICT, it is a legitimate projection like the others.

Michael

9

Glynn:

> there are many places where PROJECTION_UTM causes
> cellhd.zone to become relevant.

ah, right (eg r.info). so UTM must be a special case.
(at least for gr6)

Michael:

I don't know how this works in the underlying C code,

include/gis.h defines five fundamental projection types:

#define PROJECTION_XY 0
#define PROJECTION_UTM 1
#define PROJECTION_SP 2
#define PROJECTION_LL 3
#define PROJECTION_OTHER 99

(these can be seen in the top line of `g.region -p`)

C programs are written to behave in different ways depending on
which class of projection it is.

this all can be handled in a single method that parses
proj-parms.table except for State Plane. I can see why you
call latlon a pseudo projection, but not UTM. AFAICT, it is
a legitimate projection like the others.

Like state plane, UTM is a collection of (for want of a better
term) epsg codes. The State Plane or UTM code/zone says which
epsg-like set of parameters to use/expand.

A pure projection is a mathematical formula specifying how
geographic space bends. SP and UTM are a table of contents for
which EPSG code to use.

hope that's clearer,
Hamish

10

On Oct 16, 2009, at 10:55 PM, Hamish wrote:

A pure projection is a mathematical formula specifying how
geographic space bends. SP and UTM are a table of contents for
which EPSG code to use.

So you're referring to how GRASS handles UTM projections, not how they are actually defined?

Michael

11

From a previous post I realize r.terraflow isnt working in 6.4 yet.., so

thought I would try a dumbed-down-DEM with r.sim.water ... DEM, mask,
slope, aspect, etc look good; everything seem to run, but...no outputs
are generated.... hmmm...?
( I can take this to users' list if this is the wrong place, but seems a
lot of things are in motion here...)

thanks in advance!
Chris
----------------------------------------------------------------
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > r.sim.water -t elevin=DEM100m@PERMANENT
dxin=dx@PERMANENT dyin=dy@PERMANENT rain_val=50 infil_val=0.0
manin_val=0.1 depth=water_depth@PERMANENT disch=water_discharge@PERMANENT
err=error_map@PERMANENT niter=10 outiter=2 diffc=0.8 hmax=0.3 halpha=4.0
hbeta=0.5 --overwrite
default nwalk=68487324, rwalk=68487324.000000

Min elevation = 1.00 m
Max elevation = 3109.00 m
Mean Source Rate (rainf. excess or sediment) = 0.000002 m/s or kg/m2s
Mean flow velocity = 0.628966 m/s
Mean Mannings = 0.606325
Number of iterations = 3 cells
Time step = 39.75 s
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > g.gisenv
LOCATION_NAME=Haiti
MAPSET=PERMANENT
DIGITIZER=none
GISDBASE=/home/cgn/GRASS
MONITOR=x7
GRASS_GUI=tcltk
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > g.region -p
projection: 1 (UTM)
zone: 0
datum: wgs84
ellipsoid: a=6378137 es=0.006694379990141316
north: 2338174.90627377
south: 1879574.90627377
west: 499985.21620627
east: 1246685.21620627
nsres: 100
ewres: 100
rows: 4586
cols: 7467
cells: 34243662
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > g.list type=rast
----------------------------------------------
raster files available in mapset <PERMANENT>:
DEM100m aspect dx dyy
slope_steepness
MASK basinid dxx half_basins streams
accum curvature dxy slope tcurvature
asin_visual drainage_dir dy slope_length

----------------------------------------------
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp >

12

Please try to increase number of iterations to several hundred, run it without -t
and run it on something smaller (1000x1000) - it looks like it did not get
too far when it finished the computation, but it should have given you at least
a warning or error message (it might have run out of memory).
I need to try it out, I still run it in older versions of GRASS64,

Helena

Helena Mitasova
Associate Professor
Department of Marine, Earth
and Atmospheric Sciences
1125 Jordan Hall, Campus Box 8208
North Carolina State University
Raleigh NC 27695-8208
http://skagit.meas.ncsu.edu/~helena/

On Oct 17, 2009, at 3:49 AM, cgnicholas@alamedanet.net wrote:

From a previous post I realize r.terraflow isnt working in 6.4 yet.., so

thought I would try a dumbed-down-DEM with r.sim.water ... DEM, mask,
slope, aspect, etc look good; everything seem to run, but...no outputs
are generated.... hmmm...?
( I can take this to users' list if this is the wrong place, but seems a
lot of things are in motion here...)

thanks in advance!
Chris
----------------------------------------------------------------
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > r.sim.water -t elevin=DEM100m@PERMANENT
dxin=dx@PERMANENT dyin=dy@PERMANENT rain_val=50 infil_val=0.0
manin_val=0.1 depth=water_depth@PERMANENT disch=water_discharge@PERMANENT
err=error_map@PERMANENT niter=10 outiter=2 diffc=0.8 hmax=0.3 halpha=4.0
hbeta=0.5 --overwrite
default nwalk=68487324, rwalk=68487324.000000

Min elevation = 1.00 m
Max elevation = 3109.00 m
Mean Source Rate (rainf. excess or sediment) = 0.000002 m/s or kg/m2s
Mean flow velocity = 0.628966 m/s
Mean Mannings = 0.606325
Number of iterations = 3 cells
Time step = 39.75 s
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > g.gisenv
LOCATION_NAME=Haiti
MAPSET=PERMANENT
DIGITIZER=none
GISDBASE=/home/cgn/GRASS
MONITOR=x7
GRASS_GUI=tcltk
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > g.region -p
projection: 1 (UTM)
zone: 0
datum: wgs84
ellipsoid: a=6378137 es=0.006694379990141316
north: 2338174.90627377
south: 1879574.90627377
west: 499985.21620627
east: 1246685.21620627
nsres: 100
ewres: 100
rows: 4586
cols: 7467
cells: 34243662
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp > g.list type=rast
----------------------------------------------
raster files available in mapset <PERMANENT>:
DEM100m aspect dx dyy
slope_steepness
MASK basinid dxx half_basins streams
accum curvature dxy slope tcurvature
asin_visual drainage_dir dy slope_length

----------------------------------------------
[Raster MASK present]
GRASS 6.4.0RC5 (Haiti):~/tmp >

_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev