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

GRASS GIS GRASS GIS developer
1

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

2

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

3

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.

I'm not sure what you're using to perform the projection, but r.proj
uses the region settings which the user specifies. If you set a region
with rectangular pixels, you get rectangular pixels.

Admittedly, there isn't any automatic mechanism to adjust the region
bounds to get square pixels; you have to calculate the bounds
yourself. If you don't do this explicitly, you're likely to end up
with rectangular pixels.

[If you use "g.region res=...", GRASS will divide the region extents
by the resolution to get the number of rows and columns. These values
will then be rounded to the nearest integer, and the actual n-s and
e-w resolutions computed by dividing the extents by these integers.]

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