[GRASS-dev] ascii export and import, large file problem

GRASS GIS GRASS GIS developer
1

Helena, I think the issue when I tried this before is that grass doesn't know that x and y are in degrees while z is in meters. It either assumes they are the same or uses the zfactor option to multiply the z value by something to get it into the same units as x and y.

But if I could tell grass that the z units were meters explicitly, couldn't grass do the conversion automatically? Especially since the conversion changes as you move north.

Jerry

---- Original message ----

Date: Thu, 12 Apr 2007 21:02:37 -0400
From: Helena Mitasova <hmitaso@unity.ncsu.edu>
Subject: Re: [GRASS-dev] ascii export and import, large file problem
To: Gerald Nelson <gnelson@uiuc.edu>
Cc: grass-dev list <grass-dev@grass.itc.it>

Jerry,

r.slope.aspect should work with lat-long, it calls G_distance that says

This routine computes the distance, in meters, from
* <b>x1</b>,<b>y1</b> to <b>x2</b>,<b>y2</b>. If the projection is
* latitude-longitude, this distance is measured along the geodesic. Two
* routines perform geodesic distance calculations.

It should also support the wrap-around.

Maybe this should be added to the man page.

Helena

Helena Mitasova
Dept. of Marine, Earth and Atm. Sciences
1125 Jordan Hall, NCSU Box 8208,
Raleigh NC 27695
http://skagit.meas.ncsu.edu/~helena/

On Apr 12, 2007, at 8:10 PM, Gerald Nelson wrote:

This is related to my earlier question (to which noone responded)
about why the slope calculation can't take lat-long data and just
figure out the slope from the z values.

The srtm elevation data comes in lat-long values in cells that are
square. In order to get slope, which we need for a cost map, we
project to a utm value in the center of the region we need
(UTM37N). Grass does this projection and generates rectangular
pixels. Then we do the slope calculation and other things to
generate a cost surface.

We need the ascii output to read the cost data into our custom
neural net program (because we don't have any C programmers, just a
java programmer).

So its possible that what we observed when we did the export/import
process is a function of the square/rectangular pixel issue
interacting with the arc export/import that assumes pixels are square.

My ideal would be to not have to take the data out of lat-long, but
I have to do that to use the slope calculations.

Hope this all makes some sense.

Regards, Jerry

---- Original message ----

Date: Thu, 12 Apr 2007 23:58:16 +0100
From: Glynn Clements <glynn@gclements.plus.com>
Subject: RE: [GRASS-dev] ascii export and import, large file problem
To: "Jerry Nelson" <gnelson@uiuc.edu>
Cc: "'grass-dev'" <grass-dev@grass.itc.it>

Jerry Nelson wrote:

I'm using grass6.3 updated today so the large file support for
the ascii
commands is included. I export a file using r.out.arc and then
read it back
in using r.in.arc. The attached jpg shows the original raster
on the right.
The screen on the left is the original raster minus the
exported and
imported version. The bottom two thirds or so of the left
raster is zero, as
it should be, but the top 1/3 has a bunch of small values
(range is - to
+2.9).

My first guess is that the export->import process is changing the
vertical extent of the map slightly, so the calculation in the
upper
portion of the map is using cells which are off by one row.

What does r.info say about the bounds of the two maps?

To provide more info,

The 'after' info

Rows: 21048

Res: 119.047796

The 'before' info

Rows: 21048

Res: 119.05225396

119.05225396 - 119.047796 = 0.00445796
0.00445796 * 21048 = 93.8311421

So, the imported map has shrunk by almost a whole cell. That would
certainly explain the results.

Ah, I see where the problem lies:

The 'before' info

Res: 119.05225396
Res: 119.04779557

Your cells aren't square, but the ArcGrid format doesn't appear to
allow for non-square cells (single "cellsize" value rather than
separate x/y values). r.out.arc uses the horizontal resolution for
the
cellsize value; if the vertical resolution is different, you lose.

This specific issue can't be fixed. However, if the original data had
square cells, something is going wrong on the initial import.

We might want to add a check for this to r.out.arc. We can't actually
do anything beyond warn you that exporting will lose information,
although that's better than nothing.

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

Gerald Nelson
Professor, Dept. of Agricultural and Consumer Economics
University of Illinois, Urbana-Champaign
office: 217-333-6465
cell: 217-390-7888
315 Mumford Hall
1301 W. Gregory
Urbana, IL 61801

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

Gerald Nelson
Professor, Dept. of Agricultural and Consumer Economics
University of Illinois, Urbana-Champaign
office: 217-333-6465
cell: 217-390-7888
315 Mumford Hall
1301 W. Gregory
Urbana, IL 61801

2

On Apr 13, 2007, at 8:39 AM, Gerald Nelson wrote:

Helena, I think the issue when I tried this before is that grass doesn't know that x and y are in degrees while z is in meters. It either assumes they are the same or uses the zfactor option to multiply the z value by something to get it into the same units as x and y.

that issue relates to feet - r.slope.aspect converts x,y, or rather the horizontal distance derived from x,y
to meters (so if it is lat/long, feet, or anything else defined in PROJ_UNIT it converts it to meters)
si if you elevations are in anything else than meters - you have to be careful and convert also your elevation to meters
because GRASS has no way of knowing you elevation units. So if your elevation is in meters,
you should be fine. Let me know if not. I always project the lat/long data because I work in smaller
areas and it is easier to work with x,y than long.lat so I have never tried computation of slope
in long/lat, but by looking at the code it should work.

But if I could tell grass that the z units were meters explicitly, couldn't grass do the conversion automatically?

as I explained above - GRASS assumes it is in meters.

Especially since the conversion changes as you move north.

yes, it uses geodesic to covert - that is what I was trying to say in my email below:

"If the projection is latitude-longitude, this distance is measured along the geodesic."

I hope this helps - it is kind of messy,

Helena

Jerry

---- Original message ----

Date: Thu, 12 Apr 2007 21:02:37 -0400
From: Helena Mitasova <hmitaso@unity.ncsu.edu>
Subject: Re: [GRASS-dev] ascii export and import, large file problem
To: Gerald Nelson <gnelson@uiuc.edu>
Cc: grass-dev list <grass-dev@grass.itc.it>

Jerry,

r.slope.aspect should work with lat-long, it calls G_distance that says

This routine computes the distance, in meters, from
* <b>x1</b>,<b>y1</b> to <b>x2</b>,<b>y2</b>. Two
* routines perform geodesic distance calculations.

It should also support the wrap-around.

Maybe this should be added to the man page.

Helena

Helena Mitasova
Dept. of Marine, Earth and Atm. Sciences
1125 Jordan Hall, NCSU Box 8208,
Raleigh NC 27695
http://skagit.meas.ncsu.edu/~helena/

On Apr 12, 2007, at 8:10 PM, Gerald Nelson wrote:

This is related to my earlier question (to which noone responded)
about why the slope calculation can't take lat-long data and just
figure out the slope from the z values.

The srtm elevation data comes in lat-long values in cells that are
square. In order to get slope, which we need for a cost map, we
project to a utm value in the center of the region we need
(UTM37N). Grass does this projection and generates rectangular
pixels. Then we do the slope calculation and other things to
generate a cost surface.

We need the ascii output to read the cost data into our custom
neural net program (because we don't have any C programmers, just a
java programmer).

So its possible that what we observed when we did the export/import
process is a function of the square/rectangular pixel issue
interacting with the arc export/import that assumes pixels are square.

My ideal would be to not have to take the data out of lat-long, but
I have to do that to use the slope calculations.

Hope this all makes some sense.

Regards, Jerry

---- Original message ----

Date: Thu, 12 Apr 2007 23:58:16 +0100
From: Glynn Clements <glynn@gclements.plus.com>
Subject: RE: [GRASS-dev] ascii export and import, large file problem
To: "Jerry Nelson" <gnelson@uiuc.edu>
Cc: "'grass-dev'" <grass-dev@grass.itc.it>

Jerry Nelson wrote:

I'm using grass6.3 updated today so the large file support for
the ascii
commands is included. I export a file using r.out.arc and then
read it back
in using r.in.arc. The attached jpg shows the original raster
on the right.
The screen on the left is the original raster minus the
exported and
imported version. The bottom two thirds or so of the left
raster is zero, as
it should be, but the top 1/3 has a bunch of small values
(range is - to
+2.9).

My first guess is that the export->import process is changing the
vertical extent of the map slightly, so the calculation in the
upper
portion of the map is using cells which are off by one row.

What does r.info say about the bounds of the two maps?

To provide more info,

The 'after' info

Rows: 21048

Res: 119.047796

The 'before' info

Rows: 21048

Res: 119.05225396

119.05225396 - 119.047796 = 0.00445796
0.00445796 * 21048 = 93.8311421

So, the imported map has shrunk by almost a whole cell. That would
certainly explain the results.

Ah, I see where the problem lies:

The 'before' info

Res: 119.05225396
Res: 119.04779557

Your cells aren't square, but the ArcGrid format doesn't appear to
allow for non-square cells (single "cellsize" value rather than
separate x/y values). r.out.arc uses the horizontal resolution for
the
cellsize value; if the vertical resolution is different, you lose.

This specific issue can't be fixed. However, if the original data had
square cells, something is going wrong on the initial import.

We might want to add a check for this to r.out.arc. We can't actually
do anything beyond warn you that exporting will lose information,
although that's better than nothing.

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

Gerald Nelson
Professor, Dept. of Agricultural and Consumer Economics
University of Illinois, Urbana-Champaign
office: 217-333-6465
cell: 217-390-7888
315 Mumford Hall
1301 W. Gregory
Urbana, IL 61801

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

Gerald Nelson
Professor, Dept. of Agricultural and Consumer Economics
University of Illinois, Urbana-Champaign
office: 217-333-6465
cell: 217-390-7888
315 Mumford Hall
1301 W. Gregory
Urbana, IL 61801

3

Gerald Nelson wrote:

Helena, I think the issue when I tried this before is that grass
doesn't know that x and y are in degrees while z is in meters. It
either assumes they are the same or uses the zfactor option to
multiply the z value by something to get it into the same units as x
and y.

But if I could tell grass that the z units were meters explicitly,
couldn't grass do the conversion automatically? Especially since the
conversion changes as you move north.

http://article.gmane.org/gmane.comp.gis.grass.devel/19353/
http://article.gmane.org/gmane.comp.gis.grass.devel/19523/

--relevant quotes--
[1...]
For a real y-axis label I think we need to add a new (optional) "units"
element to the raster format (e.g. in cell_misc/).
[2...]

The user should be able to add a title for the legend.

For vlegends that's doable. For raster legends I hope to add
G_set_raster_units() and G_get_raster_units() to write/read a simple
string containing raster data units (set from "r.support units=") which
will be stored in $MAPSET/cell_misc/$RASTERMAP/units. Raster legends and
things like lat/lon NVIZ and r.shaded.relief could use the tag if it
existed.
--endquote--

apparently cell_misc/ is evil, so that moves it to
$MAPSET/units/$RASTERMAP until GRASS7 when it moves to
$MAPSET/raster/$RASTERMAP/units

I don't like to add another element dir to $MAPSET, but...

Hamish

4

Hamish wrote:

> Helena, I think the issue when I tried this before is that grass
> doesn't know that x and y are in degrees while z is in meters. It
> either assumes they are the same or uses the zfactor option to
> multiply the z value by something to get it into the same units as x
> and y.
>
> But if I could tell grass that the z units were meters explicitly,
> couldn't grass do the conversion automatically? Especially since the
> conversion changes as you move north.

http://article.gmane.org/gmane.comp.gis.grass.devel/19353/
http://article.gmane.org/gmane.comp.gis.grass.devel/19523/

--relevant quotes--
[1...]
For a real y-axis label I think we need to add a new (optional) "units"
element to the raster format (e.g. in cell_misc/).
[2...]
> The user should be able to add a title for the legend.

For vlegends that's doable. For raster legends I hope to add
G_set_raster_units() and G_get_raster_units() to write/read a simple
string containing raster data units (set from "r.support units=") which
will be stored in $MAPSET/cell_misc/$RASTERMAP/units. Raster legends and
things like lat/lon NVIZ and r.shaded.relief could use the tag if it
existed.
--endquote--

apparently cell_misc/ is evil, so that moves it to
$MAPSET/units/$RASTERMAP until GRASS7 when it moves to
$MAPSET/raster/$RASTERMAP/units

I don't like to add another element dir to $MAPSET, but...

Given that cell_misc exists and can't easily be removed right now, you
may as well use it for now.

Just use the _misc functions which I recently added, i.e.

  char *G__file_name_misc(char *, const char *, const char *, const char *, const char *);
  char *G_find_file_misc(const char *, const char *, char *, const char *);
  char *G_find_file2_misc(const char *, const char *, const char *, const char *);
  int G__make_mapset_element_misc (const char *, const char *);
  int G_open_new_misc(const char *, const char *, const char *);
  int G_open_old_misc(const char *, const char *, const char *, const char *);
  int G_open_update_misc(const char *, const char *, const char *);
  FILE *G_fopen_new_misc(const char *, const char *, const char *);
  FILE *G_fopen_old_misc(const char *, const char *, const char *, const char *);
  FILE *G_fopen_append_misc(const char *, const char *, const char *);
  FILE *G_fopen_modify_misc(const char *, const char *, const char *);
  int G_remove_misc (const char *, const char *, const char *);

rather than using sprintf(element, "cell_misc/%s", name) hacks.

As for what to do in 7.x: having one directory per map with a file for
each element is neater, but that limits the number of maps in a mapset
to the filesystem's hard link count limit (32000 for ext2/ext3). We've
had at least one person complain about this.

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

5

Hello Hamish
Thanks for bringing this up again. ISTR being too busy to reply straight away the last time and then forgetting...

On Wed, 18 Apr 2007, Hamish wrote:
[...]

--relevant quotes--
[1...]
For a real y-axis label I think we need to add a new (optional) "units"
element to the raster format (e.g. in cell_misc/).

OK this seems like a very good idea. Several times Helena has suggested the idea of adding "vertical datum" support to a location, but I was always quite negative about it in pointing out that there was not even any way of indicating that a particular raster map contained elevation data. But in fact, the obvious answer seems now to be that it should be done on a per-raster map basis rather than a per-location basis. I.e. if the metadata for a raster map indicates that it contains elevation data (i.e. units are a length measure) then there should also be a facility to specify the vertical datum the length measure is relative to.

[2...]

The user should be able to add a title for the legend.

For vlegends that's doable. For raster legends I hope to add
G_set_raster_units() and G_get_raster_units() to write/read a simple
string containing raster data units (set from "r.support units=") which
will be stored in $MAPSET/cell_misc/$RASTERMAP/units. Raster legends and
things like lat/lon NVIZ and r.shaded.relief could use the tag if it
existed.
--endquote--

I think it's worth putting enough thought in so that the feature is easily usable for more than just legend titles. Especially how to specify vertical datums. I suppose there are a few very standard ones, like WGS84 ellipsoidal height and that kind of thing. Actually I think the EPSG co-ordinate system database also includes vertical datums. Maybe it has some kind of unique codes that we could re-use. Maybe we should have a vdatum.table file with descriptions. Or maybe we should just put a human-readable name in and leave it for manual interpretation.

I'm really not too sure.

Paul

apparently cell_misc/ is evil, so that moves it to
$MAPSET/units/$RASTERMAP until GRASS7 when it moves to
$MAPSET/raster/$RASTERMAP/units

I don't like to add another element dir to $MAPSET, but...

Hamish

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

6

Paul Kelly wrote:

Hello Hamish
Thanks for bringing this up again. ISTR being too busy to reply
straight away the last time and then forgetting...

On Wed, 18 Apr 2007, Hamish wrote:
[...]
> --relevant quotes--
> [1...]
> For a real y-axis label I think we need to add a new (optional)
> "units" element to the raster format (e.g. in cell_misc/).

OK this seems like a very good idea. Several times Helena has
suggested the idea of adding "vertical datum" support to a location,
but I was always quite negative about it in pointing out that there
was not even any way of indicating that a particular raster map
contained elevation data. But in fact, the obvious answer seems now
to be that it should be done on a per-raster map basis rather than a
per-location basis. I.e. if the metadata for a raster map indicates
that it contains elevation data (i.e. units are a length measure)
then there should also be a facility to specify the vertical datum
the length measure is relative to.

> [2...]
>> The user should be able to add a title for the legend.
>
> For vlegends that's doable. For raster legends I hope to add
> G_set_raster_units() and G_get_raster_units() to write/read a simple
> string containing raster data units (set from "r.support units=")
> which will be stored in $MAPSET/cell_misc/$RASTERMAP/units. Raster
> legends and things like lat/lon NVIZ and r.shaded.relief could use
> the tag if it existed.
> --endquote--

I think it's worth putting enough thought in so that the feature is
easily usable for more than just legend titles. Especially how to
specify vertical datums. I suppose there are a few very standard
ones, like WGS84 ellipsoidal height and that kind of thing. Actually
I think the EPSG co-ordinate system database also includes vertical
datums. Maybe it has some kind of unique codes that we could re-use.
Maybe we should have a vdatum.table file with descriptions. Or maybe
we should just put a human-readable name in and leave it for manual
interpretation.

I'm really not too sure.

ok, proposal:

add raster map specific cell_misc/units and cell_misc/vertical_datum
files.

Each would hold a single line string containing the info.
You could set the value with e.g. "r.support units="
You could query with e.g. "r.info -u"

read/write using new cell_misc libgis fns described by Glynn in an
earlier email.

both fully optional, and non-existing by default, so fully backwards
compatible with old maps.

for use by modules, just start with strcmp(tolower(*string), "meters")
for hinting. units can be anything, so I prefer to allow freeform
metadata there, not just something from a list. I fear that vertical
datums may use very local names so while a common table would be helpful
we need to allow a custom entries. I'm not too sure about them though.

Hamish

7

Hamish,

can you please add your suggestion for adding cell_misc/units to
an appropriate place (development) on wiki so that it does not get burried
in the mail archives?

I am finding the unit info essential especially for modeling (I started to add
it to the file names for lack of other options) - for example, precipitation
map can be mm/month or mm/hr, inputs or outputs such as infiltration rate,
erodibility, critical sheer stress, all can have different units,
you can have temperature in deg Celsius or Fahrenheit
- there are many examples. Also the modules should be able
to check the unit (e.g. r.slope.aspect for feet or meters) and/or at least inform the user
in the --help what units are expected (I am just adding it to simwe - we forgot
to describe it in manual and there was no way to find out that
e.g. rainfall was expected in m/s rather than more common mm/hr).

Thanks, Helena

On Apr 20, 2007, at 3:28 AM, Hamish wrote:

Paul Kelly wrote:

Hello Hamish
Thanks for bringing this up again. ISTR being too busy to reply
straight away the last time and then forgetting...

On Wed, 18 Apr 2007, Hamish wrote:
[...]

--relevant quotes--
[1...]
For a real y-axis label I think we need to add a new (optional)
"units" element to the raster format (e.g. in cell_misc/).

OK this seems like a very good idea. Several times Helena has
suggested the idea of adding "vertical datum" support to a location,
but I was always quite negative about it in pointing out that there
was not even any way of indicating that a particular raster map
contained elevation data. But in fact, the obvious answer seems now
to be that it should be done on a per-raster map basis rather than a
per-location basis. I.e. if the metadata for a raster map indicates
that it contains elevation data (i.e. units are a length measure)
then there should also be a facility to specify the vertical datum
the length measure is relative to.

[2...]

The user should be able to add a title for the legend.

For vlegends that's doable. For raster legends I hope to add
G_set_raster_units() and G_get_raster_units() to write/read a simple
string containing raster data units (set from "r.support units=")
which will be stored in $MAPSET/cell_misc/$RASTERMAP/units. Raster
legends and things like lat/lon NVIZ and r.shaded.relief could use
the tag if it existed.
--endquote--

I think it's worth putting enough thought in so that the feature is
easily usable for more than just legend titles. Especially how to
specify vertical datums. I suppose there are a few very standard
ones, like WGS84 ellipsoidal height and that kind of thing. Actually
I think the EPSG co-ordinate system database also includes vertical
datums. Maybe it has some kind of unique codes that we could re-use.
Maybe we should have a vdatum.table file with descriptions. Or maybe
we should just put a human-readable name in and leave it for manual
interpretation.

I'm really not too sure.

ok, proposal:

add raster map specific cell_misc/units and cell_misc/vertical_datum
files.

Each would hold a single line string containing the info.
You could set the value with e.g. "r.support units="
You could query with e.g. "r.info -u"

read/write using new cell_misc libgis fns described by Glynn in an
earlier email.

both fully optional, and non-existing by default, so fully backwards
compatible with old maps.

for use by modules, just start with strcmp(tolower(*string), "meters")
for hinting. units can be anything, so I prefer to allow freeform
metadata there, not just something from a list. I fear that vertical
datums may use very local names so while a common table would be helpful
we need to allow a custom entries. I'm not too sure about them though.

Hamish