[GRASS-dev] r.in.gdal problems

GRASS GRASS developer
,
1

Hi,

I imported a large SRTM grid with r.in.gdal; when starting grass64 with this location there is an error.
Some part of Grass seems to work despite this error, except d.rast shows a blank screen. The east and west (see below) have the same value, except E/W. Could this cause troubles?

GRASS 6.4.0RC4 (northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america > Error in startup script: divide by zero
     while executing
"expr round($parts(w)/$parts(ewres))*$parts(ewres)"
     (procedure "MapCanvas::zoom_gregion" line 17)
     invoked from within
"MapCanvas::zoom_gregion $mon"
     (procedure "MapCanvas::create" line 40)
     invoked from within
"MapCanvas::create"
     (procedure "Gm::startmon" line 11)
     invoked from within
"Gm::startmon"
     (procedure "Gm::create" line 79)
     invoked from within
"Gm::create"
     (procedure "main" line 30)
     invoked from within
"main $argc $argv"
     (file "/usr/local/grass/grass-6.4.0RC4/etc/gm/gm.tcl" line 566)

GRASS 6.4.0RC4 (northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america > g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 61:00:01.5N
south: 9:59:58.5N
west: 179:59:58.5E
east: 179:59:58.5W
nsres: 0:00:03
ewres: 0:00:00.000007
rows: 61201
cols: 432001
cells: 26438893201

GRASS 6.4.0RC4 (northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america > r.info na_mosaic
  +----------------------------------------------------------------------------+
  | Layer: na_mosaic Date: Tue Jun 9 21:54:39 2009 |
  | Mapset: PERMANENT Login of Creator: ltoma |
  | Location: northamericaSRTM3 |
  | DataBase: /mnt/research/gis/DATA/grassdb/ |
  | Title: ( na_mosaic ) |
  | Timestamp: none |
  |----------------------------------------------------------------------------|
  | |
  | Type of Map: raster Number of Categories: 5573 |
  | Data Type: CELL |
  | Rows: 61201 |
  | Columns: 432001 |
  | Total Cells: 26438893201 |
  | Projection: Latitude-Longitude |
  | N: 61:00:01.5N S: 9:59:58.5N Res: 0:00:03 |
  | E: 179:59:58.5W W: 179:59:58.5E Res: 0:00:00.000007 |
  | Range of data: min = -291 max = 5573 |
  | |
  | Data Description: |
  | generated by r.in.gdal |
  | |
  | Comments: |
  | r.in.gdal input="na_mosaic.vrt" output="na_mosaic" |
  | |
  +----------------------------------------------------------------------------+

-Laura

2

Laura Toma wrote:

Hi,

I imported a large SRTM grid with r.in.gdal; when starting grass64
with this location there is an error.
Some part of Grass seems to work despite this error, except d.rast
shows a blank screen. The east and west (see below) have the same
value, except E/W. Could this cause troubles?

Both values are smaller than 180, so they should not give problems in a
ll projection. What could give problems is the ew resolution of
0:00:00.000007 which is an amazing resolution for a DEM covering the
whole world, especially SRTM (native resolution is 0:00:03 for the world
coverage dataset). Should that not be the same like the ns resolution,
0:00:03? Maybe something went wrong with r.in.gdal? Does gdalinfo report
the same ew resolution? The resolution of 0:00:00.000007 can cause the
"divide by zero" error when trying to set the current region and/or
display parameters. A workaround may be to start grass in text mode and
set the region resolution to 0:00:03, both ns and ew, then start the
gui. Or manually edit the WIND and DEFAULT_WIND files before launching
grass.

Markus M

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america
> Error in startup script: divide by zero
    while executing
"expr round($parts(w)/$parts(ewres))*$parts(ewres)"
    (procedure "MapCanvas::zoom_gregion" line 17)
    invoked from within
"MapCanvas::zoom_gregion $mon"
    (procedure "MapCanvas::create" line 40)
    invoked from within
"MapCanvas::create"
    (procedure "Gm::startmon" line 11)
    invoked from within
"Gm::startmon"
    (procedure "Gm::create" line 79)
    invoked from within
"Gm::create"
    (procedure "main" line 30)
    invoked from within
"main $argc $argv"
    (file "/usr/local/grass/grass-6.4.0RC4/etc/gm/gm.tcl" line 566)

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america
> g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 61:00:01.5N
south: 9:59:58.5N
west: 179:59:58.5E
east: 179:59:58.5W
nsres: 0:00:03
ewres: 0:00:00.000007
rows: 61201
cols: 432001
cells: 26438893201

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america
> r.info na_mosaic
+----------------------------------------------------------------------------+

| Layer: na_mosaic Date: Tue Jun 9 21:54:39
2009 |
| Mapset: PERMANENT Login of Creator:
ltoma |
| Location:
northamericaSRTM3 |
| DataBase:
/mnt/research/gis/DATA/grassdb/ |
| Title: ( na_mosaic
) |
| Timestamp:
none |
|----------------------------------------------------------------------------|

|
|
| Type of Map: raster Number of Categories:
5573 |
| Data Type:
CELL |
| Rows:
61201 |
| Columns:
432001 |
| Total Cells:
26438893201 |
| Projection:
Latitude-Longitude |
| N: 61:00:01.5N S: 9:59:58.5N Res:
0:00:03 |
| E: 179:59:58.5W W: 179:59:58.5E Res:
0:00:00.000007 |
| Range of data: min = -291 max =
5573 |
|
|
| Data
Description: |
| generated by
r.in.gdal |
|
|
|
Comments:
|
| r.in.gdal input="na_mosaic.vrt"
output="na_mosaic" |
|
|
+----------------------------------------------------------------------------+

-Laura_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

3

Laura Toma wrote:

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america
> g.region -p
projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 61:00:01.5N
south: 9:59:58.5N
west: 179:59:58.5E
east: 179:59:58.5W
nsres: 0:00:03
ewres: 0:00:00.000007
rows: 61201
cols: 432001
cells: 26438893201

The EW extends and number of columns would give an ewres of 0:00:03 not
0:00:00.000007. Something went wrong with the region setting and ewres
calculation of the SRTM raster, see also your r.info output below. In
this case, hack also the corresponding file in cellhd? Not really a
solution though.

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america
> r.info na_mosaic
+----------------------------------------------------------------------------+

| Layer: na_mosaic Date: Tue Jun 9 21:54:39
2009 |
| Mapset: PERMANENT Login of Creator:
ltoma |
| Location:
northamericaSRTM3 |
| DataBase:
/mnt/research/gis/DATA/grassdb/ |
| Title: ( na_mosaic
) |
| Timestamp:
none |
|----------------------------------------------------------------------------|

|
|
| Type of Map: raster Number of Categories:
5573 |
| Data Type:
CELL |
| Rows:
61201 |
| Columns:
432001 |
| Total Cells:
26438893201 |
| Projection:
Latitude-Longitude |
| N: 61:00:01.5N S: 9:59:58.5N Res:
0:00:03 |
| E: 179:59:58.5W W: 179:59:58.5E Res:
0:00:00.000007 |
| Range of data: min = -291 max =
5573 |
|
|
| Data
Description: |
| generated by
r.in.gdal |
|
|
|
Comments:
|
| r.in.gdal input="na_mosaic.vrt"
output="na_mosaic" |
|
|
+----------------------------------------------------------------------------+

-Laura_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

4

Markus,

Here is what the cellhd contains (it seems correct).

proj: 3
zone: 0
north: 61:00:01.5N
south: 9:59:58.5N
east: 179:59:58.5W
west: 179:59:58.5E
cols: 432001
rows: 61201
e-w resol: 0:00:03
n-s resol: 0:00:03
format: 3
compressed: 1

  but r.info gives:
    Projection: Latitude-Longitude |
  | N: 61:00:01.5N S: 9:59:58.5N Res: 0:00:03 |
  | E: 179:59:58.5W W: 179:59:58.5E Res: 0:00:00.000007 |

so for some reason the ewres is messed up.

Here is the tif file that I used as input in r.in.gdal:

gdalinfo na_mosaic.tif
Driver: GTiff/GeoTIFF
Files: na_mosaic.tif
Size is 432001, 61201
Coordinate System is:
GEOGCS["WGS 84",
     DATUM["WGS_1984",
         SPHEROID["WGS 84",6378137,298.2572235630016,
             AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0],
     UNIT["degree",0.0174532925199433],
     AUTHORITY["EPSG","4326"]]
Origin = (-180.000416666666666,61.000416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
   AREA_OR_POINT=Area
Image Structure Metadata:
   INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-180.0004167, 61.0004167) (180d 0'1.50"W, 61d 0'1.50"N)
Lower Left (-180.0004167, 9.9995833) (180d 0'1.50"W, 9d59'58.50"N)
Upper Right ( 180.0004167, 61.0004167) (180d 0'1.50"E, 61d 0'1.50"N)
Lower Right ( 180.0004167, 9.9995833) (180d 0'1.50"E, 9d59'58.50"N)
Center ( -0.0000000, 35.5000000) ( 0d 0'0.00"W, 35d30'0.00"N)
Band 1 Block=432001x1 Type=Int16, ColorInterp=Gray
   NoData Value=-32768
[ltoma@dover:~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america]$

On Jun 10, 2009, at 4:10 PM, Markus GRASS wrote:

Laura Toma wrote:

Hi,

I imported a large SRTM grid with r.in.gdal; when starting grass64
with this location there is an error.
Some part of Grass seems to work despite this error, except d.rast
shows a blank screen. The east and west (see below) have the same
value, except E/W. Could this cause troubles?

Both values are smaller than 180, so they should not give problems in a
ll projection. What could give problems is the ew resolution of
0:00:00.000007 which is an amazing resolution for a DEM covering the
whole world, especially SRTM (native resolution is 0:00:03 for the world
coverage dataset). Should that not be the same like the ns resolution,
0:00:03? Maybe something went wrong with r.in.gdal? Does gdalinfo report
the same ew resolution? The resolution of 0:00:00.000007 can cause the
"divide by zero" error when trying to set the current region and/or
display parameters. A workaround may be to start grass in text mode and
set the region resolution to 0:00:03, both ns and ew, then start the
gui. Or manually edit the WIND and DEFAULT_WIND files before launching
grass.

Markus M

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america

Error in startup script: divide by zero

   while executing
"expr round($parts(w)/$parts(ewres))*$parts(ewres)"
   (procedure "MapCanvas::zoom_gregion" line 17)
   invoked from within
"MapCanvas::zoom_gregion $mon"
   (procedure "MapCanvas::create" line 40)
   invoked from within
"MapCanvas::create"
   (procedure "Gm::startmon" line 11)
   invoked from within
"Gm::startmon"
   (procedure "Gm::create" line 79)
   invoked from within
"Gm::create"
   (procedure "main" line 30)
   invoked from within
"main $argc $argv"
   (file "/usr/local/grass/grass-6.4.0RC4/etc/gm/gm.tcl" line 566)

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america

g.region -p

projection: 3 (Latitude-Longitude)
zone: 0
datum: wgs84
ellipsoid: wgs84
north: 61:00:01.5N
south: 9:59:58.5N
west: 179:59:58.5E
east: 179:59:58.5W
nsres: 0:00:03
ewres: 0:00:00.000007
rows: 61201
cols: 432001
cells: 26438893201

GRASS 6.4.0RC4
(northamericaSRTM3):~/gis/DATA/SRTM/SRTM-V2-frome0srp01u.ecs.nasa.gov/SRTM3/north_america

r.info na_mosaic

+----------------------------------------------------------------------------+

| Layer: na_mosaic Date: Tue Jun 9 21:54:39
2009 |
| Mapset: PERMANENT Login of Creator:
ltoma |
| Location:
northamericaSRTM3 |
| DataBase:
/mnt/research/gis/DATA/grassdb/ |
| Title: ( na_mosaic
) |
| Timestamp:
none |
|----------------------------------------------------------------------------|

|
| Type of Map: raster Number of Categories:
5573 |
| Data Type:
CELL |
| Rows:
61201 |
| Columns:
432001 |
| Total Cells:
26438893201 |
| Projection:
Latitude-Longitude |
| N: 61:00:01.5N S: 9:59:58.5N Res:
0:00:03 |
| E: 179:59:58.5W W: 179:59:58.5E Res:
0:00:00.000007 |
| Range of data: min = -291 max =
5573 |
|
| Data
Description: |
| generated by
r.in.gdal |
|
Comments:
|
| r.in.gdal input="na_mosaic.vrt"
output="na_mosaic" |
|
+----------------------------------------------------------------------------+

-Laura_______________________________________________
grass-dev mailing list
grass-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-dev

5

Laura Toma wrote:

Markus,

Here is what the cellhd contains (it seems correct).

proj: 3
zone: 0
north: 61:00:01.5N
south: 9:59:58.5N
east: 179:59:58.5W
west: 179:59:58.5E

That's the problem. g.region or some underlying function thinks that
this raster starts at 179:59:58.5E and extends towards east until
179:59:58.5W, giving a total extend of 0:00:03 degrees longitude instead
of 363 degrees longitude. Longitude values larger than 180E or smaller
than 180E are somewhere wrapped around. As Hamish said, the safest would
probably be "gdaltranslate -srcwin 1 0 431999 61201" to chop off the
first and last column. Not so sure if it is sufficient to chop off only
the first or only the last column. Then r.in.gdal again.

BTW, the current east and west bounds are perfectly legal in grass and
allow a raster to smoothly cross the 180E/W border. So you could create
another DEM with e.g. east=90:00:01.5W and west=90:00:01.5E going over
the 180E/W border.

Good luck with that DEM:-)

Markus M

cols: 432001
rows: 61201
e-w resol: 0:00:03
n-s resol: 0:00:03
format: 3
compressed: 1

but r.info gives:
   Projection: Latitude-Longitude |
| N: 61:00:01.5N S: 9:59:58.5N Res:
0:00:03 |
| E: 179:59:58.5W W: 179:59:58.5E Res:
0:00:00.000007 |

so for some reason the ewres is messed up.