[GRASS-user] Imported .tif maps lose all data

rshepard · September 25, 2021, 2:43pm

I have 77 LiDAR topo files in .tif format. What gdalinfo tells me about a
representative file:
$ gdalinfo columbia_2010_e_dtm_35.tif
Driver: GTiff/GeoTIFF
Files: columbia_2010_e_dtm_35.tif
Size is 10745, 15264
Coordinate System is:
PROJCRS["NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet",
     BASEGEOGCRS["GCS_North_American_1983_HARN",
         DATUM["North_American_1983_HARN",
             ELLIPSOID["GRS_1980",6378137,298.257222101,
                 LENGTHUNIT["metre",1,
                     ID["EPSG",9001]]]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433,
                 ID["EPSG",9122]]]],
     CONVERSION["Lambert Conic Conformal (2SP)",
         METHOD["Lambert Conic Conformal (2SP)",
             ID["EPSG",9802]],
         PARAMETER["Latitude of false origin",45.3333333333333,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8821]],
         PARAMETER["Longitude of false origin",-120.5,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8822]],
         PARAMETER["Latitude of 1st standard parallel",45.8333333333333,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8823]],
         PARAMETER["Latitude of 2nd standard parallel",47.3333333333333,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8824]],
         PARAMETER["Easting at false origin",500000,
             LENGTHUNIT["metre",1],
             ID["EPSG",8826]],
         PARAMETER["Northing at false origin",0,
             LENGTHUNIT["metre",1],
             ID["EPSG",8827]]],
     CS[Cartesian,2],
         AXIS["easting",east,
             ORDER[1],
             LENGTHUNIT["US survey foot",0.304800609601219,
                 ID["EPSG",9003]]],
         AXIS["northing",north,
             ORDER[2],
             LENGTHUNIT["US survey foot",0.304800609601219,
                 ID["EPSG",9003]]]]
Data axis to CRS axis mapping: 1,2
Origin = (1512164.000000000000000,152311.000000000000000)
Pixel Size = (3.000000000000000,-3.000000000000000)
Metadata:
   AREA_OR_POINT=Area
   DataType=Generic
Image Structure Metadata:
   COMPRESSION=LZW
   INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 1512164.000, 152311.000) (121d 0' 8.57"W, 45d44'59.58"N)
Lower Left ( 1512164.000, 106519.000) (121d 0' 4.46"W, 45d37'27.52"N)
Upper Right ( 1544399.000, 152311.000) (120d52'34.01"W, 45d45' 1.35"N)
Lower Right ( 1544399.000, 106519.000) (120d52'30.94"W, 45d37'29.29"N)
Center ( 1528281.500, 129415.000) (120d56'19.49"W, 45d41'14.50"N)
Band 1 Block=512x256 Type=Float32, ColorInterp=Gray
   Min=159.316 Max=1037.520
   Minimum=159.316, Maximum=1037.520, Mean=262.041, StdDev=157.494
   NoData Value=-3.40282306073709653e+38
   Overviews: 5373x7632, 2687x3816, 1344x1908, 672x954, 336x477, 168x239
   Metadata:
     BandName=Band_1
     RepresentationType=ATHEMATIC
     STATISTICS_COVARIANCES=24804.31297971113
     STATISTICS_MAXIMUM=1037.5201416016
     STATISTICS_MEAN=262.04051816209
     STATISTICS_MINIMUM=159.31565856934
     STATISTICS_SKIPFACTORX=1
     STATISTICS_SKIPFACTORY=1
     STATISTICS_STDDEV=157.49385060919

I started GRASS with a new location:
grass -c /path/to/project/data/topography/columbia_2010_e_dtm_35.tif /data/grassdata/columbia_river_dtms/topography

and imported all raw data; e.g.,

r.in.gdal
in=/path/to/project/data/topography/columbia_2010_e_dtm_35.tif out=dtm_4 --o

The problem is the imported maps have no elevation data. For example,
r.report dtm_35
100%
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT |
|LOCATION: new_nevins_dock Sat Sep 25 07:28:25 2021|
|-----------------------------------------------------------------------------|
| north: 106729 east: 1544180 |
|REGION south: 60934 west: 1511876 |
| res: 3 res: 3 |
|-----------------------------------------------------------------------------|
|MASK: none |
|-----------------------------------------------------------------------------|
|MAP: (untitled) (dtm_35 in topography) |
|-----------------------------------------------------------------------------|
| Category Information | cell| % |
|#|description | count| cover|
|-----------------------------------------------------------------------------|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . . .|164373520|100.00|
|-----------------------------------------------------------------------------|
|TOTAL |164373520|100.00|
+-----------------------------------------------------------------------------+

r.stats -l in=dtm_35 out=-
100%
* no data

The data file sizes are all in the low 100Ks. I could upload an example to
fileconvoy.com if that's helpful.

What have I done incorrectly that the imported maps have no elevation data
while the raw .tif files do?

TIA,

Rich

Micha_Silver1 · September 25, 2021, 3:22pm

Rich:

Just to check: did you set the current region to the imported raster before running r.report?

From the man page:

Note that the statistics is always computed in the current geographical region.

From the Y coordinates in the gdalinfo output and the region settings reported by r.report, it seems that there is almost no overlap. Maybe that's the reason the r.report is not showing anything?

On 9/25/21 5:43 PM, Rich Shepard wrote:

I have 77 LiDAR topo files in .tif format. What gdalinfo tells me about a
representative file:
$ gdalinfo columbia_2010_e_dtm_35.tif

Corner Coordinates:
Upper Left ( 1512164.000, 152311.000) (121d 0' 8.57"W, 45d44'59.58"N)
Lower Left ( 1512164.000, 106519.000) (121d 0' 4.46"W, 45d37'27.52"N)
Upper Right ( 1544399.000, 152311.000) (120d52'34.01"W, 45d45' 1.35"N)
Lower Right ( 1544399.000, 106519.000) (120d52'30.94"W, 45d37'29.29"N)

r.in.gdal
in=/path/to/project/data/topography/columbia_2010_e_dtm_35.tif out=dtm_4 --o

Also, in your example, you imported to a GRASS raster named dtm_4, but you then run r.report on a raster named dtm_35

??

The problem is the imported maps have no elevation data. For example,
r.report dtm_35
100%
+-----------------------------------------------------------------------------+

| RASTER MAP CATEGORY REPORT |
|LOCATION: new_nevins_dock Sat Sep 25 07:28:25 2021|
|-----------------------------------------------------------------------------|

| north: 106729 east: 1544180 |
|REGION south: 60934 west: 1511876 |
| res: 3 res: 3 |
|-----------------------------------------------------------------------------|

|MASK: none |
|-----------------------------------------------------------------------------|

|MAP: (untitled) (dtm_35 in topography) |
|-----------------------------------------------------------------------------|

| Category Information | cell| % |
|#|description | count| cover|
|-----------------------------------------------------------------------------|

|*|no data. . . . . . . . . . . . . . . . . . . . . . . . . .|164373520|100.00|
|-----------------------------------------------------------------------------|

|TOTAL |164373520|100.00|
+-----------------------------------------------------------------------------+

r.stats -l in=dtm_35 out=-
100%
* no data

r.stats also honors the current region

The data file sizes are all in the low 100Ks. I could upload an example to
fileconvoy.com if that's helpful.

What have I done incorrectly that the imported maps have no elevation data
while the raw .tif files do?

TIA,

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

--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918

Carlos_Guano_Grohman · September 25, 2021, 4:07pm

How is your region settings? Do they match the DEMs? Maybe importing with r.import will help, as it can override region settings

Carlos

On Sat, 25 Sep 2021, 11:43 Rich Shepard, <rshepard@appl-ecosys.com> wrote:

I have 77 LiDAR topo files in .tif format. What gdalinfo tells me about a
representative file:
$ gdalinfo columbia_2010_e_dtm_35.tif
Driver: GTiff/GeoTIFF
Files: columbia_2010_e_dtm_35.tif
Size is 10745, 15264
Coordinate System is:
PROJCRS[“NAD_1983_HARN_StatePlane_Washington_South_FIPS_4602_Feet”,
BASEGEOGCRS[“GCS_North_American_1983_HARN”,
DATUM[“North_American_1983_HARN”,
ELLIPSOID[“GRS_1980”,6378137,298.257222101,
LENGTHUNIT[“metre”,1,
ID[“EPSG”,9001]]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433,
ID[“EPSG”,9122]]]],
CONVERSION[“Lambert Conic Conformal (2SP)”,
METHOD[“Lambert Conic Conformal (2SP)”,
ID[“EPSG”,9802]],
PARAMETER[“Latitude of false origin”,45.3333333333333,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8821]],
PARAMETER[“Longitude of false origin”,-120.5,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8822]],
PARAMETER[“Latitude of 1st standard parallel”,45.8333333333333,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8823]],
PARAMETER[“Latitude of 2nd standard parallel”,47.3333333333333,
ANGLEUNIT[“degree”,0.0174532925199433],
ID[“EPSG”,8824]],
PARAMETER[“Easting at false origin”,500000,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8826]],
PARAMETER[“Northing at false origin”,0,
LENGTHUNIT[“metre”,1],
ID[“EPSG”,8827]]],
CS[Cartesian,2],
AXIS[“easting”,east,
ORDER[1],
LENGTHUNIT[“US survey foot”,0.304800609601219,
ID[“EPSG”,9003]]],
AXIS[“northing”,north,
ORDER[2],
LENGTHUNIT[“US survey foot”,0.304800609601219,
ID[“EPSG”,9003]]]]
Data axis to CRS axis mapping: 1,2
Origin = (1512164.000000000000000,152311.000000000000000)
Pixel Size = (3.000000000000000,-3.000000000000000)
Metadata:
AREA_OR_POINT=Area
DataType=Generic
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 1512164.000, 152311.000) (121d 0’ 8.57"W, 45d44’59.58"N)
Lower Left ( 1512164.000, 106519.000) (121d 0’ 4.46"W, 45d37’27.52"N)
Upper Right ( 1544399.000, 152311.000) (120d52’34.01"W, 45d45’ 1.35"N)
Lower Right ( 1544399.000, 106519.000) (120d52’30.94"W, 45d37’29.29"N)
Center ( 1528281.500, 129415.000) (120d56’19.49"W, 45d41’14.50"N)
Band 1 Block=512x256 Type=Float32, ColorInterp=Gray
Min=159.316 Max=1037.520
Minimum=159.316, Maximum=1037.520, Mean=262.041, StdDev=157.494
NoData Value=-3.40282306073709653e+38
Overviews: 5373x7632, 2687x3816, 1344x1908, 672x954, 336x477, 168x239
Metadata:
BandName=Band_1
RepresentationType=ATHEMATIC
STATISTICS_COVARIANCES=24804.31297971113
STATISTICS_MAXIMUM=1037.5201416016
STATISTICS_MEAN=262.04051816209
STATISTICS_MINIMUM=159.31565856934
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=157.49385060919

I started GRASS with a new location:
grass -c /path/to/project/data/topography/columbia_2010_e_dtm_35.tif /data/grassdata/columbia_river_dtms/topography

and imported all raw data; e.g.,

r.in.gdal
in=/path/to/project/data/topography/columbia_2010_e_dtm_35.tif out=dtm_4 --o

The problem is the imported maps have no elevation data. For example,
r.report dtm_35
100%
±----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT |

LOCATION: new_nevins_dock Sat Sep 25 07:28:25 2021

north: 106729 east: 1544180

REGION south: 60934 west: 1511876

res: 3 res: 3

-----------------------------------------------------------------------------

MASK: none

-----------------------------------------------------------------------------

MAP: (untitled) (dtm_35 in topography)

-----------------------------------------------------------------------------

Category Information

#

-----------------------------------------------------------------------------

*

-----------------------------------------------------------------------------

TOTAL

±----------------------------------------------------------------------------+

r.stats -l in=dtm_35 out=-
100%

no data

The data file sizes are all in the low 100Ks. I could upload an example to
fileconvoy.com if that’s helpful.

What have I done incorrectly that the imported maps have no elevation data
while the raw .tif files do?

TIA,

Rich

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

rshepard · September 25, 2021, 4:41pm

On Sat, 25 Sep 2021, Micha Silver wrote:

Just to check: did you set the current region to the imported raster
before running r.report?

Micha,

Just did set the region to dtm_35 and re-ran r.report:
+-----------------------------------------------------------------------------+
| RASTER MAP CATEGORY REPORT |
|LOCATION: new_nevins_dock Sat Sep 25 09:35:06 2021|
|-----------------------------------------------------------------------------|
| north: 106729 east: 1544180 |
|REGION south: 60934 west: 1511876 |
| res: 3 res: 3 |
|-----------------------------------------------------------------------------|
|MASK: none |
|-----------------------------------------------------------------------------|
|MAP: (untitled) (dtm_35 in topography) |
|-----------------------------------------------------------------------------|
| Category Information | cell| % |
|#|description | count| cover|
|-----------------------------------------------------------------------------|
|*|no data. . . . . . . . . . . . . . . . . . . . . . . . . .|164373520|100.00|
|-----------------------------------------------------------------------------|
|TOTAL |164373520|100.00|
+-----------------------------------------------------------------------------+

and r.stats:

r.stats -l in=dtm_35 out=stats.txt
* no data

The data's still AWOL.

Also, in your example, you imported to a GRASS raster named dtm_4, but you then run r.report on a raster named dtm_35

That's becaue I copied from the input script and forgot to change it to
dtm_35 which is what I used for the examples.

Thanks,

Rich

rshepard · September 25, 2021, 4:43pm

On Sat, 25 Sep 2021, Carlos Henrique Grohmann de Carvalho wrote:

How is your region settings? Do they match the DEMs? Maybe importing with
r.import will help, as it can override region settings

Carlos,

Initially it was the default region. But, I changed the region to a test
map, dtm_35, and re-ran r.report and r.stats. There's still no data in any
of the cells.

Thanks,

Rich

Micha_Silver1 · September 25, 2021, 6:08pm

···

On 9/25/21 7:41 PM, Rich Shepard wrote:

On Sat, 25 Sep 2021, Micha Silver wrote:

Just to check: did you set the current region to the imported raster
before running r.report?

Micha,

Just did set the region to dtm_35 and re-ran r.report:

The data’s still AWOL.

Can we go back to your suggestion to post the original data somewhere?

Thanks,

Rich

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

-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918

Carlos_Guano_Grohman · September 25, 2021, 6:16pm

Rich

I believe r.in.gdal will require the region to match each of the tifs you’re importing. Try using r.import with the option extent set to “input”, so it won’t do region cropping

C

···

Prof. Carlos Henrique Grohmann
Institute of Energy and Environment - Univ. of São Paulo, Brazil

Digital Terrain Analysis | GIS | Remote Sensing -

http://carlosgrohmann.com
http://orcid.org/0000-0001-5073-5572

Can’t stop the signal.

rshepard · September 25, 2021, 6:29pm

On Sat, 25 Sep 2021, Carlos Henrique Grohmann de Carvalho wrote:

I believe r.in.gdal will require the region to match each of the tifs
you're importing. Try using r.import with the option extent set to
"input", so it won't do region cropping

Carlos,

All .tif files in that directory have the same region. They cover the
Columbia River from the mouth to Bonneville Dam. But, those 146 miles are
separated into 77 separate reach files.

Regards,

Rich

Micha_Silver1 · September 25, 2021, 6:32pm

Hi Carlos:

···

On 9/25/21 9:16 PM, Carlos Henrique Grohmann de Carvalho wrote:

Rich

I believe r.in.gdal will require the region to match each of the tifs you’re importing. Try using r.import with the option extent set to

“input”, so it won’t do region cropping

I don’t think that is correct. r.in.gdal is one of the (very few) raster modules that does NOT honor the current region settings. You can import a raster from “anywhere”, but then, as always, you must set the current region to match that raster for further analysis.

What’s more, in your previous post you mentioned that r.import “overrides region settings”. IMHO this is also not correct. r.import will reproject an input dataset to current Location’s coordinate system, if necessary. Regarding the region settings, by default with the extent parameter left at “input”, the whole input raster will be imported. If you set the the extent parameter to “region” then only that subset of the input raster that falls in the current region will be imported. But in both cases, the current region is not altered.

Best,

Micha

C

On Sat, Sep 25, 2021 at 1:43 PM Rich Shepard <rshepard@appl-ecosys.com> wrote:

On Sat, 25 Sep 2021, Carlos Henrique Grohmann de Carvalho wrote:

How is your region settings? Do they match the DEMs? Maybe importing with
r.import will help, as it can override region settings

Carlos,

Initially it was the default region. But, I changed the region to a test
map, dtm_35, and re-ran r.report and r.stats. There’s still no data in any
of the cells.

Thanks,

Rich

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

–

Prof. Carlos Henrique Grohmann
Institute of Energy and Environment - Univ. of São Paulo, Brazil

Digital Terrain Analysis | GIS | Remote Sensing -

http://carlosgrohmann.com
http://orcid.org/0000-0001-5073-5572

Can’t stop the signal.
_______________________________________________
grass-user mailing list
[grass-user@lists.osgeo.org](mailto:grass-user@lists.osgeo.org)
[https://lists.osgeo.org/mailman/listinfo/grass-user](https://lists.osgeo.org/mailman/listinfo/grass-user)

-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918

Carlos_Guano_Grohman · September 25, 2021, 8:44pm

Hi Micha

You’re right. That’s what happens when I try to answer quickly while keeping an eye on a 5-year old. Thanks.

Carlos

···

Prof. Carlos Henrique Grohmann
Institute of Energy and Environment - Univ. of São Paulo, Brazil

Digital Terrain Analysis | GIS | Remote Sensing -

http://carlosgrohmann.com
http://orcid.org/0000-0001-5073-5572

Can’t stop the signal.