[GRASSLIST:9481] SRTM shifted 0.5 the res in X and Y axis from the r.mapcalc output

GRASS GRASS user
1

This might be connected with my [GRASSLIST:9409].

SRTM tiles imported with r.in.srtm are shifted 0.5 the resolution in X
and Y axis from the native Grass rasters cretaed with r.mapcalc in the
same location.

See the attached picture and notes below. I'm not doing anything
unusuall. Why are these two rasters shifted? I don't think they should
be. I guess that when importing data like SRTM (BIL), where the
coordinates are of the cell centre, into Grass, where coordinates are of
the cell ll corner, shuch data should shifted to conform to Grass grid.
Correct me please if this is rubbish what I suppose, but I'm puzzled
that I have two rasters of the same res in Grass and they don't fit each
other.

# 1. import a 3" SRTM tile
r.in.srtm input=/home/maciek/gis/dane/srtm2/N49E023.hgt.zip
output=N49E023

g.region rast=N49E023 -p projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:01.5N
south: 48:59:58.5N
west: 22:59:58.5E
east: 24:00:01.5E
nsres: 0:00:03
ewres: 0:00:03
rows: 1201
cols: 1201

# colour it and zoom to the top-right corner of the tile to see single
# cells
r.colors N49E023 col=grey.eq
d.rast N49E023
d.zoom

# create a GRASS raster there, in the same resolution
r.mapcalc 'test=rand(0,1000)'

# increase the resolution so that any small shift is visible
g.region res=0:0:0.1 -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:09N
south: 49:59:48N
west: 23:59:54E
east: 24:00:12E
nsres: 0:00:00.1
ewres: 0:00:00.1
rows: 210
cols: 180

# display the Grass raster and overlay an SRTM
d.erase; d.rast test; d.rast N49E023 -o

# here's the shift

See the attachment - grey is the SRTM, r.mapcalc ouput in random colors.

Maciek

--------------------
,Kuchnia polska PWE" - Idealny prezent pod choinkê!
Najs³ynniejsza ksi±¿ka kucharska, która pomo¿e na nowo odkryæ przysmaki polskiej kuchni
Z okazji ¶wi±t 10% zni¿ki na stronie http://www.pwe.com.pl/51552.xml

(attachments)

SRTMvsGRASSshift.png

2

Maciek,

since we have spent significant time on GDAL vs SRTM
vs BIL sv GRASS cell positioning, I wonder that
SRTM in GRASS is shifted. It's at least unlikely.
See the r.in.srtm script for comments.

Probably you are facing a nearest-neighbor problem?
I don't trust d.zoom to much...

Markus

On Wed, Dec 14, 2005 at 08:28:51PM +0100, Maciek Sieczka wrote:

This might be connected with my [GRASSLIST:9409].

SRTM tiles imported with r.in.srtm are shifted 0.5 the resolution in X
and Y axis from the native Grass rasters cretaed with r.mapcalc in the
same location.

See the attached picture and notes below. I'm not doing anything
unusuall. Why are these two rasters shifted? I don't think they should
be. I guess that when importing data like SRTM (BIL), where the
coordinates are of the cell centre, into Grass, where coordinates are of
the cell ll corner, shuch data should shifted to conform to Grass grid.
Correct me please if this is rubbish what I suppose, but I'm puzzled
that I have two rasters of the same res in Grass and they don't fit each
other.

# 1. import a 3" SRTM tile
r.in.srtm input=/home/maciek/gis/dane/srtm2/N49E023.hgt.zip
output=N49E023

g.region rast=N49E023 -p projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:01.5N
south: 48:59:58.5N
west: 22:59:58.5E
east: 24:00:01.5E
nsres: 0:00:03
ewres: 0:00:03
rows: 1201
cols: 1201

# colour it and zoom to the top-right corner of the tile to see single
# cells
r.colors N49E023 col=grey.eq
d.rast N49E023
d.zoom

# create a GRASS raster there, in the same resolution
r.mapcalc 'test=rand(0,1000)'

# increase the resolution so that any small shift is visible
g.region res=0:0:0.1 -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:09N
south: 49:59:48N
west: 23:59:54E
east: 24:00:12E
nsres: 0:00:00.1
ewres: 0:00:00.1
rows: 210
cols: 180

# display the Grass raster and overlay an SRTM
d.erase; d.rast test; d.rast N49E023 -o

# here's the shift

See the attachment - grey is the SRTM, r.mapcalc ouput in random colors.

Maciek

--------------------
,Kuchnia polska PWE" - Idealny prezent pod choink?!
Najs?ynniejsza ksi??ka kucharska, kt?ra pomo?e na nowo odkry? przysmaki polskiej kuchni
Z okazji ?wi?t 10% zni?ki na stronie http://www.pwe.com.pl/51552.xml

--
Markus Neteler <neteler itc it> http://mpa.itc.it
ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
MPBA - Predictive Models for Biol. & Environ. Data Analysis
Via Sommarive, 18 - 38050 Povo (Trento), Italy

3

I don't think there is a problem. I worried about and checked this when I started converting and filling SRTM data.

Since the cell centers of a SRTM tile range from degree to the next degree (so there is a 1 cell overlap between all tiles), the GRASS region will appear to be off by 1.5 seconds - it won't line up on even degrees.

The problem, in your example at least may be that you used d.zoom before r.mapcalc. d.zoom doesn't restrain itself to cell boundaries (as you can see from the region after zooming), and r.mapcalc will create the test raster at whatever region is set. Then, of course they will not align.

On Dec 14, 2005, at 1:28 PM, Maciek Sieczka wrote:

This might be connected with my [GRASSLIST:9409].

SRTM tiles imported with r.in.srtm are shifted 0.5 the resolution in X
and Y axis from the native Grass rasters cretaed with r.mapcalc in the
same location.

See the attached picture and notes below. I'm not doing anything
unusuall. Why are these two rasters shifted? I don't think they should
be. I guess that when importing data like SRTM (BIL), where the
coordinates are of the cell centre, into Grass, where coordinates are of
the cell ll corner, shuch data should shifted to conform to Grass grid.
Correct me please if this is rubbish what I suppose, but I'm puzzled
that I have two rasters of the same res in Grass and they don't fit each
other.

# 1. import a 3" SRTM tile
r.in.srtm input=/home/maciek/gis/dane/srtm2/N49E023.hgt.zip
output=N49E023

g.region rast=N49E023 -p projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:01.5N
south: 48:59:58.5N
west: 22:59:58.5E
east: 24:00:01.5E
nsres: 0:00:03
ewres: 0:00:03
rows: 1201
cols: 1201

# colour it and zoom to the top-right corner of the tile to see single
# cells
r.colors N49E023 col=grey.eq
d.rast N49E023
d.zoom

# create a GRASS raster there, in the same resolution
r.mapcalc 'test=rand(0,1000)'

# increase the resolution so that any small shift is visible
g.region res=0:0:0.1 -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:09N
south: 49:59:48N
west: 23:59:54E
east: 24:00:12E
nsres: 0:00:00.1
ewres: 0:00:00.1
rows: 210
cols: 180

# display the Grass raster and overlay an SRTM
d.erase; d.rast test; d.rast N49E023 -o

# here's the shift

See the attachment - grey is the SRTM, r.mapcalc ouput in random colors.

-----
William Kyngesburye <kyngchaos@kyngchaos.com>
http://www.kyngchaos.com/

"This is a question about the past, is it? ... How can I tell that the past isn't a fiction designed to account for the discrepancy between my immediate physical sensations and my state of mind?"

- The Ruler of the Universe

4

Ack, I was partially right about the zoom thing. d.zoom restrains itself to the region resolution, but it's aligning to a whole degree, not the 1.5 sec offset that it started at. Thus it appears to be half a cell off.

On Dec 14, 2005, at 2:50 PM, William Kyngesburye wrote:

I don't think there is a problem. I worried about and checked this when I started converting and filling SRTM data.

Since the cell centers of a SRTM tile range from degree to the next degree (so there is a 1 cell overlap between all tiles), the GRASS region will appear to be off by 1.5 seconds - it won't line up on even degrees.

The problem, in your example at least may be that you used d.zoom before r.mapcalc. d.zoom doesn't restrain itself to cell boundaries (as you can see from the region after zooming), and r.mapcalc will create the test raster at whatever region is set. Then, of course they will not align.

-----
William Kyngesburye <kyngchaos@kyngchaos.com>
http://www.kyngchaos.com/

"History is an illusion caused by the passage of time, and time is an illusion caused by the passage of history."

- Hitchhiker's Guide to the Galaxy

5

This might be connected with my [GRASSLIST:9409].

SRTM tiles imported with r.in.srtm are shifted 0.5 the resolution in X
and Y axis from the native Grass rasters cretaed with r.mapcalc in the
same location.

..

# 1. import a 3" SRTM tile
r.in.srtm input=/home/maciek/gis/dane/srtm2/N49E023.hgt.zip
output=N49E023

g.region rast=N49E023 -p projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:01.5N
south: 48:59:58.5N
west: 22:59:58.5E
east: 24:00:01.5E
nsres: 0:00:03
ewres: 0:00:03
rows: 1201
cols: 1201

# colour it and zoom to the top-right corner of the tile to see single
# cells
r.colors N49E023 col=grey.eq
d.rast N49E023
d.zoom

Markus and William are correct to suspect d.zoom. What happens is that
you pick a spot on the screen with the mouse and d.zoom aligns itself to
the nearest 0:00:03 resolution click. Maybe this is a bug, maybe it is a
feature which needs a new flag to turn off. I don't know.

e.g. in UTM try setting your extent to 200075,220025 with res=50 and try
d.zoom. It will end up on xxxx50,xxx100.

# create a GRASS raster there, in the same resolution
r.mapcalc 'test=rand(0,1000)'

# increase the resolution so that any small shift is visible
g.region res=0:0:0.1 -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 50:00:09N
south: 49:59:48N
west: 23:59:54E
east: 24:00:12E
nsres: 0:00:00.1
ewres: 0:00:00.1
rows: 210
cols: 180

try changing the res before d.zoom and it should be ok.

regards,
Hamish

6

On czw, 2005-12-15 at 22:08 +1300, Hamish wrote:

> This might be connected with my [GRASSLIST:9409].
>
> SRTM tiles imported with r.in.srtm are shifted 0.5 the resolution in X
> and Y axis from the native Grass rasters cretaed with r.mapcalc in the
> same location.
..
> # 1. import a 3" SRTM tile
> r.in.srtm input=/home/maciek/gis/dane/srtm2/N49E023.hgt.zip
> output=N49E023
>
> g.region rast=N49E023 -p projection: 3 (Latitude-Longitude)
> zone: 0
> datum: wgs84
> ellipsoid: wgs84
> north: 50:00:01.5N
> south: 48:59:58.5N
> west: 22:59:58.5E
> east: 24:00:01.5E
> nsres: 0:00:03
> ewres: 0:00:03
> rows: 1201
> cols: 1201
>
> # colour it and zoom to the top-right corner of the tile to see single
> # cells
> r.colors N49E023 col=grey.eq
> d.rast N49E023
> d.zoom

Markus and William are correct to suspect d.zoom.
What happens is that
you pick a spot on the screen with the mouse and d.zoom aligns itself to
the nearest 0:00:03 resolution click.

Thanks Hamish, Markus and William. I think I get it now.

Maybe this is a bug, maybe it is a
feature which needs a new flag to turn off. I don't know.

IMHO d.zoom should behave in such way that it wouldn't cause this kind
of doubts as I had.

If d.zoom'ed my SRTM and reprojected this zoomed part only, my output
would become shifted - up to 3 arc sec horizontally and vertically. Same
if I wnated to do any r.mapcalc, r.to.vect, r.patch, else. This should
not take place - but it does. How are users supposed to know this
danger? Looks like quite a serious issue to me - data corruption threat.

e.g. in UTM try setting your extent to 200075,220025 with res=50 and try
d.zoom. It will end up on xxxx50,xxx100.

> # create a GRASS raster there, in the same resolution
> r.mapcalc 'test=rand(0,1000)'
>
> # increase the resolution so that any small shift is visible
> g.region res=0:0:0.1 -p
> projection: 3 (Latitude-Longitude)
> zone: 0
> datum: wgs84
> ellipsoid: wgs84
> north: 50:00:09N
> south: 49:59:48N
> west: 23:59:54E
> east: 24:00:12E
> nsres: 0:00:00.1
> ewres: 0:00:00.1
> rows: 210
> cols: 180

try changing the res before d.zoom and it should be ok.

I'm affraid I don't understand. Change res to what value? What for?

Cheers,
Maciek

--------------------
W polskim Internecie s± setki milionów stron. My przekazujemy Tobie tylko najlepsze z nich!
http://katalog.epf.pl/

7

> try changing the res before d.zoom and it should be ok.

I'm affraid I don't understand. Change res to what value? What for?

g.region res=0:00:00.1
d.zoom

instead of

d.zoom
g.region res=0:00:00.1

Hamish

8

On pi±, 2005-12-16 at 12:46 +1300, Hamish wrote:

> > try changing the res before d.zoom and it should be ok.
>
> I'm affraid I don't understand. Change res to what value? What for?

g.region res=0:00:00.1
d.zoom

instead of

d.zoom
g.region res=0:00:00.1

right

Anyway, what about the data corruption threat? Shall I fill a bug
report?

Maciek

--------------------
W polskim Internecie s± setki milionów stron. My przekazujemy Tobie tylko najlepsze z nich!
http://katalog.epf.pl/