[GRASS-user] Projection Units and Values Negatively Affect r.topidx Results

GRASS GRASS user
21

Rich Shepard wrote:

On Wed, 17 Feb 2010, Markus Neteler wrote:

I think that a change of the base unit of a location requires a new
location and reprojection of the data. I may be wrong, though.

Markus/Hamish:

  I'm missing something critical here.

  I defined a new location (called oregonM) with the following projection
information and units:

PROJ_INFO:

name: Lambert Conformal Conic
proj: lcc
datum: nad83
ellps: grs80
lat_1: 43
lat_2: 45.5
lat_0: 41.45
lon_0: -120.5
x_0: 400000
y_0: 0
no_defs: defined
nadgrids: WO

PROJ_UNITS:

unit: meter
units: meters
meters: 1

  Using the wxPython GUI I import the source DEM specifying hdr.adf as the
ARC/GRID source. On the options tab I check "override projection information
and use location's projection" because that file has the units in feet:

[rshepard@salmo /usr4/Oregon]$ less dem10-nw/prj.adf
Projection LAMBERT Datum NAD83 Zunits NO Units 3.2808400000 Spheroid GRS1980 Xshift 0.0000000000 Yshift 0.0000000000 Parameters
43 0 0.000 /* 1st standard parallel
45 30 0.000 /* 2nd standard parallel -120 30 0.000 /* central meridian
41 45 0.000 /* latitude of projection's origin 400000.00000 /* false easting (meters) 0.00000 /* false northing (meters)

  However, after the file is imported and I run r.info on the map I see:

Type of Map: raster Number of Categories: 5729 Data Type: CELL
Rows: 18724
Columns: 24657
Total Cells: 461677668
Projection: Lambert Conformal Conic
N: 1735231.43724681 S: 1120657.1324816 Res: 32.82281055
E: 952979.22944945 W: 143667.18968253 Res: 32.82281055

  I cannot seem to be able to convert from feet to meters during the import.

If I followed your post correctly, you should not be converting anything during the inport. You import into a Location defined with projection params exactly like the original ArcInfo file.
Then create a new, empty Location with units meters. Switch to the new location, import the region vector polygon (from v.in.region in the original Location), and set the region to that vector.
Now you should be ready to do r.proj to re-project the original (feet) raster to the new region with units=meters.

The command used was:
r.in.gdal -o input="/usr4/Oregon/dem10-nw/hdr.adf" output="demNW10"

  What have I missed here?

Rich

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

--
Micha Silver
Arava Development Co. +972-52-3665918
http://surfaces.co.il

22

On Thu, 18 Feb 2010, Micha Silver wrote:

If I followed your post correctly, you should not be converting anything during the inport. You import into a Location defined with projection params exactly like the original ArcInfo file.

Micha,

   I realized that I need to use r.proj (and v.proj) to make the conversions.

Then create a new, empty Location with units meters.

location=oregonM, mapset=PERMANENT PROJ_UNITS:
oregonM/PERMANENT/PROJ_UNITS unit: meter
units: meters
meters: 1

Switch to the new location, import the region vector polygon (from
v.in.region in the original Location), and set the region to that vector.

Source location:

GRASS 6.4.0svn (Oregon):/usr4/grassbase > v.in.region out=coverage type=area
Building topology for vector map <coverage>...
Registering primitives...
2 primitives registered
6 vertices registered
Building areas...
  100%
1 areas built
1 isles built

New destination location:

GRASS 6.4.0svn (oregonM):/usr4/grassbase > v.proj in=coverage loc=Oregon
map=PERMANENT out=coverage --o

Number of areas: 1
Number of isles: 1

GRASS 6.4.0svn (oregonM):/usr4/grassbase > g.region vect=coverage -p
projection: 99 (Lambert Conformal Conic)
zone: 0
datum: nad83
ellipsoid: grs80
north: 562245.46557333
south: 374923.21748089
west: 43789.75941523
east: 290468.06913619
nsres: 32.82324305
ewres: 32.82479171
rows: 5707
cols: 7515
cells: 42888105

Now you should be ready to do r.proj to re-project the original (feet)
raster to the new region with units=meters.

   That's the theory, but it just is not working. No matter how many times I
re-do the process the resolutions remain in International Feet, not meters.
ARRGGGHHH-H-H!

   The immediate question is: why is this not working? It's probably user
error, but I am not seeing what I'm doing incorrectly.

Toda,

Rich

23

Rich wrote:

> 2) use v.in.region

v.in.region out=coverage type=area

> 2) get resolution in meters with g.region -m (applies
only to raster map reprojection)

GRASS 6.4.0svn (oregonM):/usr4/grassbase > g.region -m
n=1735231.43724681
s=1120657.1324816
w=143667.18968253
e=952979.22944945
nsres=32.82281055
ewres=32.82281055
rows=18724
cols=24657
cells=461677668

The resolution is still in feet (or is this now seen
as 32.8 meters?

This is where I'm hung up.

try switching to grass 6.5 for a minute, and running r.proj -p
and or -g. The output of -g can be pasted directly into g.region,
and then it is an easy 'g.region -p', 'g.region res= -ap' to
clean it up. It is a real breath of fresh air after all the
quiting & restarting grass + v.in.region tricks, although that
is still a nice check to do.

check that the PROJ_UNITS in your source location is valid.

do not expect coords to be the same, they start off counting
from 0,0 somewhere far off to the south west and cover a
different number of units to get to Oregon based on their length.

Hamish

24

Hi,
by browsing the man pages I have found that r.transect should do the work. Howevere at the moment it works by telling the module east, north, azimuth and distance, while maybe it would be good to have also another couple of options like:

1 choose a vector make with sql query to tell the module where to sample elevations
2 be able to draw a cross sections and get and ascii file distance-elevations (this could be achieved by d.profile if there was such an option)

At present the only way I can get a topographic profile to load into a spreadsheet is to use the oldtcltk, query the raster (d.where) and then use the d.where output as the input data for r.profile (still in oldtcltk) which gives me an ascii distance-elevation file. Indetail:

cat d.where-output | r.profile input=dem output=elevation_points

Hope other people can contribute
Francesco

John Tate wrote:

Hey Francesco,

I'd be interested to hear what other people have to say about this.

I used d.what.rast as I had ground control points, from GPS, for the profiles with which to validate the DEM. Add columns and update them with d.what.rast and export (.csv), then I put that together with the distance info from the GPS data.

Hope that helps. But probably not if your are drawing random profiles.

John

On Thursday 18 February 2010 15:42:31 Francesco Mirabella wrote:

Hi all,
I am trying to use the r.profile module on a DEM. I have tried on both
grass64_rc5 and grass7 with wxpython gui.
In 64_rc5 there still exist the "old" option which allowed the user to
interactively select the end-points by mouse, but the function does not
appear to work as I get:
---------------------------------------------------------------------------
------ r.profile -i input=DEM output=profile
Using resolution [90.0161]
No socket to connect to for monitor <cairo7>.
---------------------------------------------------------------------------
------ I can use it only with the oldtcltk

In grass 7 the possibility to choose the end-points with mouse is not
there any more
can someone give hints??

By the way, the display profiler in the monitor works good in 64_rc5 -
not i grass7 (I guess it's the "old" d.profile") but I cannot see any
way to export the data in ascii format (distance-elevation) to load them
into a spreadsheet

more hints??
many thanks

Francesco

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

--
**********************************************
Francesco Mirabella,
Geologia Strutturale e Geofisica
Universita' di Perugia,
Dipartimento di Scienze della Terra,
Piazza Universita' 1, 06100 Perugia (Italy)
tel: ++39.(0)75.584.7948
fax: ++39.(0)75.585.2603
skype: francesco.mirabella
web: http://www.unipg.it/~mirabell/
**********************************************

25

On Thu, 18 Feb 2010, Hamish wrote:

try switching to grass 6.5 for a minute, and running r.proj -p and or -g.
The output of -g can be pasted directly into g.region, and then it is an
easy 'g.region -p', 'g.region res= -ap' to clean it up. It is a real
breath of fresh air after all the quiting & restarting grass + v.in.region
tricks, although that is still a nice check to do.

Hamish,

   Wow! After all of yesterday's anguish, 6.5's r.proj did the trick. That's
teriffic! So, should I continue this project in 6.5 or revert to the 6.4
revision 41119?

check that the PROJ_UNITS in your source location is valid.

   Yes: they're in English units and the destination units are metric.

do not expect coords to be the same, they start off counting from 0,0
somewhere far off to the south west and cover a different number of units
to get to Oregon based on their length.

   Besides this, the bounds of the source location are in feet while the
bounds of the destination location are in meters. Different numbers and
different magnitudes.

   Now to reproject all source data and re-do all the project-specific data.
Well, there goes the weekend, but at least I know the numbers and analyses
are correct if this case ever goes to trial and the plaintiff's attorney
challenges my analytical methods, results, and interpretations.

Many thanks again,

Rich

26

Rich:

Wow! After all of yesterday's anguish, 6.5's r.proj
did the trick. That's teriffic! So, should I continue this
project in 6.5 or revert to the 6.4 revision 41119?

the modules are the same, the only difference is the new -p
and -g flags in 6.5 which can save you some time by not needing
to use the v.in.region where-am-I trick.

the stable branch is recommended for production use.

but at least I know the numbers and analyses are correct

as mentioned earlier, 'd.grid -w' to overlay a lat/lon grid
can act as a great reality check if you have a printed map
of the area you trust & can compare against.

Hamish

27

Francesco Mirabella wrote:

I am trying to use the r.profile module on a DEM. I have
tried on both grass64_rc5 and grass7 with wxpython gui.

In 64_rc5 there still exist the "old" option which allowed
the user to interactively select the end-points by mouse,
but the function does not appear to work as I get:
--------------------------------------------------------------------------------
r.profile -i input=DEM output=profile
Using resolution [90.0161]
No socket to connect to for monitor <cairo7>.
---------------------------------------------------------------------------------
I can use it only with the oldtcltk

Right, the new GUIs do not use Xmon windows (started with d.mon)
and in general you can not use d.* and interactive xmon modules
with them. You have to use the built-in GUI tools instead.

In this case use the GUI profile tool, which is in the map
display window just to the right of the magnifying glass buttons.
The button looks like a little line graph.

In grass 7 the possibility to choose the end-points with
mouse is not there any more

as the interactive Xmonitors have been completely removed there.
(no more 'd.mon x0')

By the way, the display profiler in the monitor works good
in 64_rc5 - not i grass7 (I guess it's the "old" d.profile")

The wxGUI profile tool from the map display window button?

but I cannot see any way to export the data in ascii format
(distance-elevation) to load them into a spreadsheet

if you use the query raster map button in the map display window
GUI you will see the easting,northing,value text printed to
the layer manager output window. you can then Save output from
there or highlight+right-click Copy the end points to a text
file and then pass that to r.profile. Not not very smooth- file
a wish in the trac bug/wish system to save the profile values to
a text file from the gui profile tool if you like.

John Tate wrote:

I used d.what.rast as I had ground control points, from GPS,
for the profiles with which to validate the DEM. Add columns
and update them with d.what.rast and export (.csv), then I put
that together with the distance info from the GPS data.

v.what.rast might help, or r.what.

(I prefer v.rast.stats with a small buffer around each point to
smooth any noise; some scripts in addons for aiding with that)

Francesco:

by browsing the man pages I have found that r.transect should
do the work. Howevere at the moment it works by telling the
module east, north, azimuth and distance,

r.transect is just a wrapper around the r.profile module. All
it does is convert the heading,distance to x2,y2 for you.

while maybe it would be good to have also another couple of
options like:

1 choose a vector make with sql query to tell the module where
to sample elevations

v.what.rast

2 be able to draw a cross sections and get and ascii file
distance-elevations (this could be achieved by d.profile if
there was such an option)

isn't this what the r.profile output gives?

At present the only way I can get a topographic profile to
load into a spreadsheet is to use the oldtcltk, query the
raster (d.where) and then use the d.where output as the input
data for r.profile (still in oldtcltk) which gives me an ascii
distance-elevation file. Indetail:

cat d.where-output | r.profile input=dem output=elevation_points

fwiw you can simplify:
d.where | r.profile input=rastermap output=-

Hamish

28

On Sat, Feb 20, 2010 at 2:12 AM, Hamish <hamish_b@yahoo.com> wrote:

Rich:

Wow! After all of yesterday's anguish, 6.5's r.proj
did the trick. That's teriffic! So, should I continue this
project in 6.5 or revert to the 6.4 revision 41119?

the modules are the same, the only difference is the new -p
and -g flags in 6.5 which can save you some time by not needing
to use the v.in.region where-am-I trick.

Please let's backport it for 6.4.1 then.

Markus

29

Hi Hamish,
many thank you for your messege,

Hamish wrote:

Francesco Mirabella wrote:

I am trying to use the r.profile module on a DEM. I have
tried on both grass64_rc5 and grass7 with wxpython gui.

In 64_rc5 there still exist the "old" option which allowed
the user to interactively select the end-points by mouse,
but the function does not appear to work as I get:
--------------------------------------------------------------------------------
r.profile -i input=DEM output=profile
Using resolution [90.0161]
No socket to connect to for monitor <cairo7>.
---------------------------------------------------------------------------------
I can use it only with the oldtcltk

Right, the new GUIs do not use Xmon windows (started with d.mon)
and in general you can not use d.* and interactive xmon modules
with them. You have to use the built-in GUI tools instead.

Ok, this is what I rechoned

In this case use the GUI profile tool, which is in the map
display window just to the right of the magnifying glass buttons.
The button looks like a little line graph.

In grass 7 the possibility to choose the end-points with
mouse is not there any more

as the interactive Xmonitors have been completely removed there.
(no more 'd.mon x0')

By the way, the display profiler in the monitor works good
in 64_rc5 - not i grass7 (I guess it's the "old" d.profile")

The wxGUI profile tool from the map display window button?

Yes, it is really a nice tool with beatifull graphical rendering

but I cannot see any way to export the data in ascii format
(distance-elevation) to load them into a spreadsheet

if you use the query raster map button in the map display window
GUI you will see the easting,northing,value text printed to
the layer manager output window. you can then Save output from
there or highlight+right-click Copy the end points to a text
file and then pass that to r.profile. Not not very smooth- file
a wish in the trac bug/wish system to save the profile values to
a text file from the gui profile tool if you like.

John Tate wrote:

I used d.what.rast as I had ground control points, from GPS,
for the profiles with which to validate the DEM. Add columns
and update them with d.what.rast and export (.csv), then I put
that together with the distance info from the GPS data.

v.what.rast might help, or r.what.

(I prefer v.rast.stats with a small buffer around each point to
smooth any noise; some scripts in addons for aiding with that)

Francesco:

by browsing the man pages I have found that r.transect should
do the work. Howevere at the moment it works by telling the
module east, north, azimuth and distance,

r.transect is just a wrapper around the r.profile module. All
it does is convert the heading,distance to x2,y2 for you.

while maybe it would be good to have also another couple of
options like:

1 choose a vector make with sql query to tell the module where
to sample elevations

v.what.rast

2 be able to draw a cross sections and get and ascii file
distance-elevations (this could be achieved by d.profile if
there was such an option)

isn't this what the r.profile output gives?

Yes, however, maybe the possibility to make a straight ascii export of distance-elevation data from the wxGUI profile tool would be good. (d.profile did this before by giving plot.A, plot.B etc..). At present the user can export images (png, svg ect.) but not the data from the wxGUI profile tool. I think this would be good if one needs to plot more data along a topographic profile (in the same graph).
Of course at present the r.profile command-line workaround works good, but maybe a more intuitive tools would be good (in the wxGUI mood :-))

best wishes

Francesco

At present the only way I can get a topographic profile to
load into a spreadsheet is to use the oldtcltk, query the
raster (d.where) and then use the d.where output as the input
data for r.profile (still in oldtcltk) which gives me an ascii
distance-elevation file. Indetail:

cat d.where-output | r.profile input=dem output=elevation_points

fwiw you can simplify:
d.where | r.profile input=rastermap output=-

Hamish

30

On Monday 22 February 2010 09:22:55 Francesco Mirabella wrote:

Hi Hamish,
many thank you for your messege,

Hamish wrote:
> Francesco Mirabella wrote:
>> I am trying to use the r.profile module on a DEM. I have
>> tried on both grass64_rc5 and grass7 with wxpython gui.
>>
>> In 64_rc5 there still exist the "old" option which allowed
>> the user to interactively select the end-points by mouse,
>> but the function does not appear to work as I get:
>> ------------------------------------------------------------------------
>>-------- r.profile -i input=DEM output=profile
>> Using resolution [90.0161]
>> No socket to connect to for monitor <cairo7>.
>> ------------------------------------------------------------------------
>>--------- I can use it only with the oldtcltk
>
> Right, the new GUIs do not use Xmon windows (started with d.mon)
> and in general you can not use d.* and interactive xmon modules
> with them. You have to use the built-in GUI tools instead.

Ok, this is what I rechoned

> In this case use the GUI profile tool, which is in the map
> display window just to the right of the magnifying glass buttons.
> The button looks like a little line graph.
>
>> In grass 7 the possibility to choose the end-points with
>> mouse is not there any more
>
> as the interactive Xmonitors have been completely removed there.
> (no more 'd.mon x0')
>
>> By the way, the display profiler in the monitor works good
>> in 64_rc5 - not i grass7 (I guess it's the "old" d.profile")
>
> The wxGUI profile tool from the map display window button?

Yes, it is really a nice tool with beatifull graphical rendering

>> but I cannot see any way to export the data in ascii format
>> (distance-elevation) to load them into a spreadsheet
>
> if you use the query raster map button in the map display window
> GUI you will see the easting,northing,value text printed to
> the layer manager output window. you can then Save output from
> there or highlight+right-click Copy the end points to a text
> file and then pass that to r.profile. Not not very smooth- file
> a wish in the trac bug/wish system to save the profile values to
> a text file from the gui profile tool if you like.
>
> John Tate wrote:
>> I used d.what.rast as I had ground control points, from GPS,
>> for the profiles with which to validate the DEM. Add columns
>> and update them with d.what.rast and export (.csv), then I put
>> that together with the distance info from the GPS data.
>
> v.what.rast might help, or r.what.
>
> (I prefer v.rast.stats with a small buffer around each point to
> smooth any noise; some scripts in addons for aiding with that)

Cheers for pointing out my obvious mistake, v.what.rast is what I meant.

I looked at v.rast.stats. Nice that it updates the table but for this purpose,
it produces a lot of extra unneeded columns (shame you can't select which
stats to add to the table - or perhaps you can?).

> Francesco:
>> by browsing the man pages I have found that r.transect should
>> do the work. Howevere at the moment it works by telling the
>> module east, north, azimuth and distance,
>
> r.transect is just a wrapper around the r.profile module. All
> it does is convert the heading,distance to x2,y2 for you.
>
>> while maybe it would be good to have also another couple of
>> options like:
>>
>> 1 choose a vector make with sql query to tell the module where
>> to sample elevations
>
> v.what.rast
>
>> 2 be able to draw a cross sections and get and ascii file
>> distance-elevations (this could be achieved by d.profile if
>> there was such an option)
>
> isn't this what the r.profile output gives?

Yes, however, maybe the possibility to make a straight ascii export of
distance-elevation data from the wxGUI profile tool would be good.
(d.profile did this before by giving plot.A, plot.B etc..). At present
the user can export images (png, svg ect.) but not the data from the
wxGUI profile tool. I think this would be good if one needs to plot more
data along a topographic profile (in the same graph).
Of course at present the r.profile command-line workaround works good,
but maybe a more intuitive tools would be good (in the wxGUI mood :-))

best wishes

Francesco

>> At present the only way I can get a topographic profile to
>> load into a spreadsheet is to use the oldtcltk, query the
>> raster (d.where) and then use the d.where output as the input
>> data for r.profile (still in oldtcltk) which gives me an ascii
>> distance-elevation file. Indetail:
>>
>> cat d.where-output | r.profile input=dem output=elevation_points
>
> fwiw you can simplify:
> d.where | r.profile input=rastermap output=-
>
>
>
> Hamish

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

31

John Tate wrote:

I looked at v.rast.stats. Nice that it updates the table
but for this purpose,
it produces a lot of extra unneeded columns (shame you
can't select which
stats to add to the table - or perhaps you can?).

I don't think you can. (beyond v.db.dropcol after)
file a wish in the trac system if you like.

Hamish