[GRASSLIST:1215] Combining raster layers with r.mapcalc

GRASS GRASS user
1

Hey there:

In running through the SEEDS tutorial, one of the tasks is combining two raster layers into a single raster using r.mapcalc. The individual layers show the data properly, but the resulting combined layer contains no data. Here is what I do:

Mapset <wbdickinson> in Location <leics2>
GRASS 5.0.2 > r.mapcalc
Enter expressions, "end" when done.
mapcalc> buffers = mwaybuf + railbuf
mapcalc> end

  100%

Mapset <wbdickinson> in Location <leics2>
GRASS 5.0.2 > d.rast buffers
  100%

Mapset <wbdickinson> in Location <leics2>
GRASS 5.0.2 > r.report buffers
r.stats: 100%
+-----------------------------------------------------------------------------
| RASTER MAP CATEGORY REPORT |LOCATION: leics2 Fri Sep 12 13:30:07 2003
|-----------------------------------------------------------------------------
| north: 322000 east: 456000 |REGION south: 310000 west: 444000 | res: 50 res: 50 |-----------------------------------------------------------------------------
|MASK:none |-----------------------------------------------------------------------------
|MAP: (untitled) (buffers in wbdickinson) |-----------------------------------------------------------------------------
| Category Information |#|description |-----------------------------------------------------------------------------
|*|no data +-----------------------------------------------------------------------------

Any thoughts would be appreciated. TIA
   Bill

--

Bill Dickinson
GIS Specialist
NASA Goddard Space Flight Center
Environmental & Safety Branch, Code 250
wdickins@pop700.gsfc.nasa.gov

2

Dear Grass Users

Recently, I bought a Landsat7 ETM+ data from Michigan State University,

which in HDF format:

( l71146029_02920010925_b10.l1g, ., l71146029_02920010925_b80.l1g).

And to import this file into GRASS, I learned that I must have GDAL
installed in my PC.

So following to [Building GDAL From Source] instruction, I installed it into
my Linux system.

Following is the directory:

/usr/local/grass5

/home/gdal-1.1.9

And the instruction said that:

[In order to run GDAL after installing it is necessary for the shared
library to be findable. This can often be accomplished by setting
LD_LIBRARY_PATH to include /usr/local/lib].

But I don't know how to do this.

Please tell me how I could do this step.

Heartfully thanks to any advice.

Ghalip.Yasin (Alipu.Yasen)

Tsukuba University, Japan

s01k154@shako.sk.tsukuba.ac.jp

3

Dear Grass Users

I have a scanned map.
And I want to digitize this map onscreen use GRASS.
Could anyone tell me, please, if this is possibel on Grass, and how could I
do it?

Heartfully Thanks for any advice.

Sincerely
*********************************
Ghalip Yasin
University of Tsukuba, Japan
e-mail: s01k154@shako.sk.tsukuba.ac.jp
*********************************

4

On Fri, 12 Sep 2003, Bill Dickinson Jr wrote:

Hey there:

In running through the SEEDS tutorial, one of the tasks is combining
two raster layers into a single raster using r.mapcalc. The
individual layers show the data properly, but the resulting combined
layer contains no data. Here is what I do:

This is *probably* because the tutorial was written for GRASS 4.x where
there was no concept of a NULL value, only zero. In GRASS 5.x r.mapcalc,
where any part of the input maps contains no data, the output will have no
data. See here for a fuller discussion (and why the tutorial cannot be
updated):
http://intevation.de/rt/webrt?serial_num=1338

Maybe use r.patch instead?

Paul

5

On Fri, 12 Sep 2003, Paul Kelly wrote:

On Fri, 12 Sep 2003, Bill Dickinson Jr wrote:

> Hey there:
>
> In running through the SEEDS tutorial, one of the tasks is combining
> two raster layers into a single raster using r.mapcalc. The
> individual layers show the data properly, but the resulting combined
> layer contains no data. Here is what I do:

This is *probably* because the tutorial was written for GRASS 4.x where
there was no concept of a NULL value, only zero. In GRASS 5.x r.mapcalc,
where any part of the input maps contains no data, the output will have no
data. See here for a fuller discussion (and why the tutorial cannot be
updated):
http://intevation.de/rt/webrt?serial_num=1338

Yes, it's because the zeros in the input rasters are being converted to
NULL. Another way round is:

r.mapcalc "try1 = isnull(rail)"
r.mapcalc "rail1 = if(try1, 0, 1, 0)"

or equivalently:

r.mapcalc "rail1 = if(isnull(rail), 0, 1, 0)"

and the same for mway.

Maybe use r.patch instead?

r.patch is more elegant, though:

r.patch input=railbuf,mwaybuf output=buffers

Roger

Paul

--
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand@nhh.no

6

Bill Dickinson Jr wrote:

In running through the SEEDS tutorial, one of the tasks is combining
two raster layers into a single raster using r.mapcalc. The
individual layers show the data properly, but the resulting combined
layer contains no data. Here is what I do:

Mapset <wbdickinson> in Location <leics2>
GRASS 5.0.2 > r.mapcalc
Enter expressions, "end" when done.
mapcalc> buffers = mwaybuf + railbuf
mapcalc> end

The SEEDS tutorial was written for GRASS 4.x. One of the main changes
between 4.x and 5.x was the introduction of a distinguished "null"
(aka "no-data") value; 4.x used to use zero for this purpose (as well
as to represent actual zero values).

Whereas the 4.x version of r.buffer set cells outside of the buffer
zone to zero, the 5.x version sets them to null. The behaviour of
r.mapcalc is that, in most cases (and this is true of addition),
whenever one of the operands is null, the result is null.

So, in 4.x the above expression would result in computations such as:

  x + 0 => x
  0 + y => y
  x + y => x + y

while in 5.x it would result in:

  x + null => null
  null + y => null
  x + y => x + y

Consequently, your result map is (almost) entirely null.

Possible solutions:

1. Use r.null to first change the nulls into zeroes:

  r.null map=mwaybuf null=0
  r.null map=railbuf null=0

then use r.mapcalc as above.

2. Re-write the r.mapcalc expression as:

GRASS 5.0.2 > r.mapcalc
mapcalc> buffers = if(isnull(mwaybuf),0,mwaybuf) + if(isnull(railbuf),0,railbuf)
mapcalc> end

3. Use r.patch instead of r.mapcalc:

  r.patch input=mwaybuf,railbuf output=buffers

This last option won't produce quite the same result, although it
might actually be more useful.

--
Glynn Clements <glynn.clements@virgin.net>

7

Dear Ghalip,

Yes, you can digitize your map using GRASS. The exact steps should be
like this:

For unrectified map.
1. import the map into GRASS (ie. "r.in.tiff")
2. digitizing (v.digit)

For rectified map:
1. produce XY location within your $GRASSDATA folder (using the GRASS
start-up dialog box). Its extent is slightly higher than the dimesnions
of your map.
2. produce the TARGET location with the cart. projection equal to your
map ()using the GRASS start-up dialog box
3. In XY location rectify the map (i.group, i.target, i. points,
i.rectify)
4. Under the TARGET location open your map.
5. Digitize it (v.digit)

For details see pages 63-70:

@book{nete02,
author={Markus Neteler and Helena Mitá¹ová},
title={{Open Source GIS: A GRASS GIS Approach}},
publisher={{Kluver Academic Publishers}},
note={1st.edition},
year={2002}}

Regards,

Rado

Dòa Pi, 2003-09-12 at 14:57, Alip.Yasin napísal:

Dear Grass Users

I have a scanned map.
And I want to digitize this map onscreen use GRASS.
Could anyone tell me, please, if this is possibel on Grass, and how could I
do it?

Heartfully Thanks for any advice.

Sincerely
*********************************
Ghalip Yasin
University of Tsukuba, Japan
e-mail: s01k154@shako.sk.tsukuba.ac.jp
*********************************

--
Radoslav Bonk M.S.
Dept. of Physical Geography and Geoecology
Faculty of Sciences, Comenius University
Mlynska Dolina 842 15, Bratislava, SLOVAKIA
tel: +421 905 968 127 e-mail: rbonk@host.sk

8

Dear Rado Bonk and Jaime Carrera

Thank you for all of your advice.
Your answers very, very pleased me.
Because as I learned from you, by using GRASS in my PC, there is no any need
for me to walk a long way to go to other students Lab and use their
Arc/info.
I will try to digitize my map using v.digit under support of GRASS tutorial.
Thank you a lot.

Sincerely
Ghalip Yasin
Tsukuba University, Japan

9

And the instruction said that:

[In order to run GDAL after installing it is necessary for the shared
library to be findable. This can often be accomplished by setting
LD_LIBRARY_PATH to include /usr/local/lib].

But I don't know how to do this.

edit /etc/ld.so.conf and add /usr/local/lib to the list if is isn't
there. Then run /sbin/ldconfig.

or maybe LD_LIBRARY_PATH=/usr/local/lib:${LD_LIBRARY_PATH}
  ??

Hamish

10

Hello Hamish, and GRASS Users

Warmly thank you all for your kindly help to me.

I installed GDAL and also add /usr/local/lib to /etc/ld.so.conf followed
Mr. Hamish's advice.

> [In order to run GDAL after installing it is necessary for the shared
> library to be findable. This can often be accomplished by setting
> LD_LIBRARY_PATH to include /usr/local/lib].
edit /etc/ld.so.conf and add /usr/local/lib to the list if is isn't
there. Then run /sbin/ldconfig.

And after that I tried to import Landsat7 ETM+
l71146029_02920010925_b10.l1g,
using r.in.gdal command.
And it's said that [ Error 4, not recognized as a supported file ].
But I learned that GDAL can handle with HDF.
So I wonder why I could't import Landsat7 etm+ into GRASS.
Is my data which named l71146029_02920010925_b10.l1g was not a HDF,
because, I thought, the extensions was '"l1g" not "hdf".
How could I inroduce the "xxx...xxx...xxx.l1g" into GRASS?
(my installed GDAL verstion was gdal1.1.9)

Thank you very much for your help.

Ghalip Yasin
Tsukuba Univeristy, Japan
s01k154@shako.sk.tsukuba.ac.jp

11

On Sat, Sep 13, 2003 at 05:09:24PM +0900, Alip.Yasin wrote:

And after that I tried to import Landsat7 ETM+
l71146029_02920010925_b10.l1g,
using r.in.gdal command.
And it's said that [ Error 4, not recognized as a supported file ].
But I learned that GDAL can handle with HDF.
So I wonder why I could't import Landsat7 etm+ into GRASS.
Is my data which named l71146029_02920010925_b10.l1g was not a HDF,
because, I thought, the extensions was '"l1g" not "hdf".
How could I inroduce the "xxx...xxx...xxx.l1g" into GRASS?
(my installed GDAL verstion was gdal1.1.9)

Your Lansat is probably in FAST format. To read FAST dataset using GDAL
you should give to r.in.gdal name of the header file. Read here for more
information:

http://www.remotesensing.org/gdal/frmt_fast.html

             Andrey

--
Andrey V. Kiselev
Home phone: +7 812 5274898 ICQ# 26871517

12

Dear GRASS Users

Firstly Thanks to Mr. Andrey Kiselev and Mr. Markus Neteler for their
support to me.

I had problem with importing a Landsat7 file into GRASS.
   my PC System: Linux; GRASS 5.0.2 (April 2003); GDAL 1.1.9;
   importing file: Landsat7 ETM+; Format: EOS-HDF;
When I type command:
  GRASS:/>gdalinfo L7114029_02920010925_hdf.L1g
   Driver: HDF4/Hierarchical Data Format Release 4
   Size is 512,512
   Coordinate System is
   Subdatasets:
          SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:
"L71146029_02920010925_hdf.L1g":0
          SUBDATASET_1_DESC=[7261x8161] [L71146029-02920010925.B10 (8_bit
unsigned integer)]
            ........
   Corner Coordinates:
     Upper Left (0.0, 0.0)
    .............
     Center (256.0,256.0)

And when I type command:
  GRASS:/>gdal -config --formats
there was a lot of other data formats, but not include format HDF4.
And when I tried to import file L71146029_02920010925_hdf.L1g use command:
   GRASS:/>r.in.gdal input=L71146029_02920010925_hdf.L1g output=test1
it said:
   Segmentation Fault.
Please, tell me, what cause this problem, and how could I import my Landsat7
HDF4 data into GRASS?

Heartfully Thanks

Ghalip Yasin
Tsukuba University, Japan
s01k154@shako.sk.tsukuba.ac.jp

On Sat, Sep 13, 2003 at 05:09:24PM +0900, Alip.Yasin wrote:
> And after that I tried to import Landsat7 ETM+
> l71146029_02920010925_b10.l1g,
> using r.in.gdal command.
> And it's said that [ Error 4, not recognized as a supported file ].
> But I learned that GDAL can handle with HDF.
> So I wonder why I could't import Landsat7 etm+ into GRASS.
> Is my data which named l71146029_02920010925_b10.l1g was not a HDF,
> because, I thought, the extensions was '"l1g" not "hdf".
> How could I inroduce the "xxx...xxx...xxx.l1g" into GRASS?
> (my installed GDAL verstion was gdal1.1.9)

Your Lansat is probably in FAST format. To read FAST dataset using GDAL
you should give to r.in.gdal name of the header file. Read here for more
information:

13

On Tue, Sep 16, 2003 at 03:04:04PM +0900, Alip.Yasin wrote:

I had problem with importing a Landsat7 file into GRASS.
   my PC System: Linux; GRASS 5.0.2 (April 2003); GDAL 1.1.9;
   importing file: Landsat7 ETM+; Format: EOS-HDF;

Alip,

The following assertions are controversal:

When I type command:
  GRASS:/>gdalinfo L7114029_02920010925_hdf.L1g
   Driver: HDF4/Hierarchical Data Format Release 4
   Size is 512,512
   Coordinate System is
   Subdatasets:
          SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:"L71146029_02920010925_hdf.L1g":0
          SUBDATASET_1_DESC=[7261x8161] [L71146029-02920010925.B10 (8_bit unsigned integer)]

and

And when I type command:
  GRASS:/>gdal -config --formats
there was a lot of other data formats, but not include format HDF4.

If you don't have hdf4 compiled in, you will not see report for HDF4
dataset.

Note, that you should use

$ r.in.gdal input=HDF4_SDS:UNKNOWN:"L71146029_02920010925_hdf.L1g":0

command to import the first band into GRASS. Read HDF4 help page for
details:

  http://www.remotesensing.org/gdal/frmt_hdf4.html

            Andrey

--
Andrey V. Kiselev
Home phone: +7 812 5274898 ICQ# 26871517

14

Dear GRASS Users
Dear Andrey Kiselev
Dear Nishida Kenrou

Thank you very much for your advice.
For import my Landsat7 data I tried command:
     >gdalinfo myfilename_hdf.l1g
and
     >gdalinfo nyfilename_mtl.l1g
it goes pretty well and I got a lot of information which contained in these
files.
But when it comes to commands( I tried in X,Y & Lon.Lan coordinates, and got
the same error message):
     >r.in.gdal input=myfilename_hdf.l1g output=name1
     >r.in.gdal input=HDF4_SDS:UNKNOWN:"myfilename_hdf.l1g":0 output=name1
all it said were:
       Segmentaion Fault;
to commands:
     >r.in.bin input=myfilename_hdf.l1g output=name1 bytes=1 north=5048400
south=4830600 east=747000 west=502200 r=217800 c=244800
it said:
          Bytes do not match Filesize.
          Filesize 65535 ... Total Bytes 1777832448
          Try bytes=0 or adjusting input parameteres.
So, after this I tried every combination of parameteres which I could
thought, but all were the same result.
Please, any advice about how could I go out of this wood.

Heartfully thanks
Ghalip Yasin (Alipu Yasen)
Tsukuba University, Japan
s01k154@shako.sk.tsukuba.ac.jp

On Tue, Sep 16, 2003 at 03:04:04PM +0900, Alip.Yasin wrote:
> I had problem with importing a Landsat7 file into GRASS.
> my PC System: Linux; GRASS 5.0.2 (April 2003); GDAL 1.1.9;
> importing file: Landsat7 ETM+; Format: EOS-HDF;

Alip,

The following assertions are controversal:

> When I type command:
> GRASS:/>gdalinfo L7114029_02920010925_hdf.L1g
> Driver: HDF4/Hierarchical Data Format Release 4
> Size is 512,512
> Coordinate System is
> Subdatasets:
>

SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:"L71146029_02920010925_hdf.L1g":0

> SUBDATASET_1_DESC=[7261x8161] [L71146029-02920010925.B10

(8_bit unsigned integer)]

and

> And when I type command:
> GRASS:/>gdal -config --formats
> there was a lot of other data formats, but not include format HDF4.

If you don't have hdf4 compiled in, you will not see report for HDF4
dataset.

Note, that you should use

$ r.in.gdal input=HDF4_SDS:UNKNOWN:"L71146029_02920010925_hdf.L1g":0

command to import the first band into GRASS. Read HDF4 help page for
details:

  http://www.remotesensing.org/gdal/frmt_hdf4.html

  Andrey

--
Andrey V. Kiselev
Home phone: +7 812 5274898 ICQ# 26871517

15

Dear Dr. Kenlou NISHIDA,

Heartfully thanks to you for spent your precious time on my problem and gave
me a solution.
Followed your method, I successfully introduced my Landsat7 ETM+ HDF4 file
into GRASS.
For more once, thank you very much.

And also THANKS A LOT to:
   Mr. Andrey Kiselev, Mr. Hamish, Mr. Markus Neteler for yours kindliness
advice.

Sincerely
Ghalip Yasin (Alipu Yasen)
Tsukuba University, Japan
s01k154@shako.sk.tsukuba.ac.jp

Open a mapset with XY-coordinate.

GRASS> r.in.bin input=landsat7filename output=outputname c=xxxx r=xxxx

west=xxxx east=xxxx

south=xxxx north=xxxx byte=x