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

[hellik]

Am I
missing an option?

nothing obvious.

try the same steps as I've done (with adjusted paths) and _post the command and the command output_ here.

- grass78/grass80 -e -c columbia_2010_e_dtm_35.tif \grassdata\rastertest
- enter into the location just created
- g.region -p
- r.in.gdal input=columbia_2010_e_dtm_35.tif output=columbia_2010_e_dtm_35
- g.region -p -a raster=columbia_2010_e_dtm_35 align=columbia_2010_e_dtm_35
- r.stats -l input=columbia_2010_e_dtm_35
- r.report map=columbia_2010_e_dtm_35

Helmut

[rshepard]

On Sat, 25 Sep 2021, Helmut Kudrnovsky wrote:

try the same steps as I've done (with adjusted paths) and _post the command and the command output_ here.

- grass78/grass80 -e -c columbia_2010_e_dtm_35.tif \grassdata\rastertest
- enter into the location just created

Helmut,

Problem at this step:

$ grass -e -c columbia_2010_e_dtm_35.tif ./rastertest
Starting GRASS GIS...
Creating new GRASS GIS location <rastertest>...
Cleaning up temporary files...

# PROBLEM (PWD is /data/grassdata/):
$ grass rastertest Starting GRASS GIS...
ERROR: </data/grassdata> is not a valid GRASS Location because PERMANENT Mapset is missing
Maybe you meant a different directory.
Exiting...

$ grass /data/grassdata/rastertest/
Starting GRASS GIS...
ERROR: </data/grassdata> is not a valid GRASS Location because PERMANENT Mapset is missing
Maybe you meant a different directory.
Exiting...

# The rastertest location exists with the PERMANENT mapset:
ls rastertest/PERMANENT/
DEFAULT_WIND MYNAME PROJ_INFO PROJ_UNITS PROJ_WKT VAR WIND sqlite/

So, grass here is seeing the GISDBASE as a location, even when the location
is specified on the command line with either that mapset name or the full
path?

I've been having this issue the past few weeks.

Rich

[mlennert]

Le 25 septembre 2021 22:09:57 GMT+02:00, Rich Shepard <rshepard@appl-ecosys.com> a écrit :

On Sat, 25 Sep 2021, Helmut Kudrnovsky wrote:

try the same steps as I've done (with adjusted paths) and _post the command and the command output_ here.

- grass78/grass80 -e -c columbia_2010_e_dtm_35.tif \grassdata\rastertest
- enter into the location just created

Helmut,

Problem at this step:

$ grass -e -c columbia_2010_e_dtm_35.tif ./rastertest
Starting GRASS GIS...
Creating new GRASS GIS location <rastertest>...
Cleaning up temporary files...
Cleaning up temporary files...

# PROBLEM (PWD is /data/grassdata/):
$ grass rastertest
Starting GRASS GIS...
ERROR: </data/grassdata> is not a valid GRASS Location because PERMANENT Mapset is missing
Maybe you meant a different directory.
Exiting...

$ grass /data/grassdata/rastertest/
Starting GRASS GIS...
ERROR: </data/grassdata> is not a valid GRASS Location because PERMANENT Mapset is missing
Maybe you meant a different directory.
Exiting...

# The rastertest location exists with the PERMANENT mapset:
ls rastertest/PERMANENT/
DEFAULT_WIND MYNAME PROJ_INFO PROJ_UNITS PROJ_WKT VAR WIND sqlite/

So, grass here is seeing the GISDBASE as a location, even when the location
is specified on the command line with either that mapset name or the full
path?

I've been having this issue the past few weeks.

The startup script expects you to provide the full path to a mapset, so when you give it 'data/grassdata/rastertest' it assumes that 'rastertest' is the mapset, so 'grassdata' must be the location. It then checks for a PERMANENT mapset in that location but cannot find it.

After creating your location it only contains one single mapset: 'PERMANENT', so to start GRASS GIS, you have to use

grass /data/grassdata/rastertest/PERMANENT

Moritz

[rshepard]

On Sat, 25 Sep 2021, Helmut Kudrnovsky wrote:

try the same steps as I've done (with adjusted paths) and _post the
command and the command output_ here.
- g.region -p
- r.in.gdal input=columbia_2010_e_dtm_35.tif output=columbia_2010_e_dtm_35
- g.region -p -a raster=columbia_2010_e_dtm_35 align=columbia_2010_e_dtm_35
- r.stats -l input=columbia_2010_e_dtm_35
- r.report map=columbia_2010_e_dtm_35

Helmut,

Results in attached 574 line file. There are some cells with data, but most
without data (82%).

I can display the data with cells; for this one map they appear to cover
only the river and not the terrain north of it. But, this is only one of 77
files.

Since this test worked I need to figure out why it didn't work when I
imported all 77 .tif files using the bash script. Ah! Thought: when
r.in.gdal imports all 77 files where does GRASS get the values for setting
g.region -p?

If g.region -p displays the region of last map imported then the other 76
maps will be outside that region. Is there a module that uses all maps in
that mapset to set the region, or is that supposed to be done by default?

Thanks very much,

Rich

(attachments)

grass-raster-test.txt (30.9 KB)

[rshepard]

On Sat, 25 Sep 2021, Moritz Lennert wrote:

The startup script expects you to provide the full path to a mapset, so
when you give it 'data/grassdata/rastertest' it assumes that 'rastertest'
is the mapset, so 'grassdata' must be the location. It then checks for a
PERMANENT mapset in that location but cannot find it.

After creating your location it only contains one single mapset:
'PERMANENT', so to start GRASS GIS, you have to use
grass /data/grassdata/rastertest/PERMANENT

Moritz,

Thank you. I've forgotten that and set the PWD to the location when starting
GRASS.

I just posted the rest of the raster test and worked in the PERMANENT
mapset.

Best regards,

Rich

[mlennert]

Le 26 septembre 2021 00:50:13 GMT+02:00, Rich Shepard <rshepard@appl-ecosys.com> a écrit :

On Sat, 25 Sep 2021, Helmut Kudrnovsky wrote:

try the same steps as I've done (with adjusted paths) and _post the
command and the command output_ here.
- g.region -p
- r.in.gdal input=columbia_2010_e_dtm_35.tif output=columbia_2010_e_dtm_35
- g.region -p -a raster=columbia_2010_e_dtm_35 align=columbia_2010_e_dtm_35
- r.stats -l input=columbia_2010_e_dtm_35
- r.report map=columbia_2010_e_dtm_35

Helmut,

Results in attached 574 line file. There are some cells with data, but most
without data (82%).

I can display the data with cells; for this one map they appear to cover
only the river and not the terrain north of it. But, this is only one of 77
files.

Since this test worked I need to figure out why it didn't work when I
imported all 77 .tif files using the bash script. Ah! Thought: when
r.in.gdal imports all 77 files where does GRASS get the values for setting
g.region -p?

You feed the information to g.region with the 'raster' parameter.

-p just means pretty print the region settings, it doesn't set the region.

If g.region -p displays the region of last map imported then the other 76
maps will be outside that region. Is there a module that uses all maps in
that mapset to set the region, or is that supposed to be done by default?

If you want to set the region to the extension covered by multiple raster maps, just feed all the maps' names to the 'raster' parameter.

But in your case, while still testing, you should set the region to one map, then test with r.report, then set the region to the next map, test, etc.

Another hint: if you just want to quickly test if the imported map has any values, use r.info which is region independent. With the -r flage you get the min-max range.

Moritz

[rshepard]

On Sat, 25 Sep 2021, Moritz Lennert wrote:

You feed the information to g.region with the 'raster' parameter.
-p just means pretty print the region settings, it doesn't set the region.

Mortitz,

I know this.

If you want to set the region to the extension covered by multiple raster
maps, just feed all the maps' names to the 'raster' parameter.

Ah! I hadn't thought of this.

But in your case, while still testing, you should set the region to one
map, then test with r.report, then set the region to the next map, test,
etc.

I don't think I need to test any more, but I'll probably do one more anyway.

Another hint: if you just want to quickly test if the imported map has any
values, use r.info which is region independent. With the -r flage you get
the min-max range.

r.info is what lead me to the issue. I didn't save the output from all the
maps. I seem to recall that as the script ran there were issues of the map
being imported not being in the current region's extent.

All my previous work with GRASS since 4.0 had one large map that covered a
larger region than needed for the project so I would define a project area
within it. These LiDAR topography maps cover a very large spatial area and I
have no idea with one (or two) actually include the project's area. So I
need to import them all then find the appropriate one or two.

Thanks very much,

Rich

[Micha_Silver1]

Hi Rich:

···

On 9/26/21 2:36 AM, Rich Shepard wrote:

All my previous work with GRASS since 4.0 had one large map that covered a
larger region than needed for the project so I would define a project area
within it. These LiDAR topography maps cover a very large spatial area and I
have no idea with one (or two) actually include the project’s area. So I
need to import them all then find the appropriate one or two.

If that’s the case, perhaps a quick gdalinfo loop over all 77 topography maps in advance will help to find the relevant ones. You get the corner coordinates in the gdalinfo output, so you’ll be able to find those that overlap the project region before importing to GRASS.

Micha

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

[rshepard]

On Sun, 26 Sep 2021, Micha Silver wrote:

If that's the case, perhaps a quick `gdalinfo` loop over all 77 topography
maps in advance will help to find the relevant ones. You get the corner
coordinates in the gdalinfo output, so you'll be able to find those that
overlap the project region before importing to GRASS.

Micha,

That's an excellent idea. I was going to use r.info or r.report on all maps
and build a text table of the corners, then find the minimum and maximum of
each corner and build the region using those values.

I have the lon-lat of the project area and can reproject that to State Plane
then find the one or two maps that enclose that area.

It's common for many of us who are closely focused on an issue to not see
other possibilities. Outsiders see things we don't.

Toda,

Rich

[mlennert]

If you set your region correctly before importing all the files, you could also use the -r flag of r.in.gdal, no ?

Might possibly be slower, though, then determining which files to import before running r.in.gdal as Micha suggests

Moritz

Le 26 septembre 2021 14:44:46 GMT+02:00, Rich Shepard <rshepard@appl-ecosys.com> a écrit :

On Sun, 26 Sep 2021, Micha Silver wrote:

If that's the case, perhaps a quick `gdalinfo` loop over all 77 topography
maps in advance will help to find the relevant ones. You get the corner
coordinates in the gdalinfo output, so you'll be able to find those that
overlap the project region before importing to GRASS.

Micha,

That's an excellent idea. I was going to use r.info or r.report on all maps
and build a text table of the corners, then find the minimum and maximum of
each corner and build the region using those values.

I have the lon-lat of the project area and can reproject that to State Plane
then find the one or two maps that enclose that area.

It's common for many of us who are closely focused on an issue to not see
other possibilities. Outsiders see things we don't.

Toda,

Rich

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

[rshepard]

On Sun, 26 Sep 2021, Moritz Lennert wrote:

If you set your region correctly before importing all the files, you could
also use the -r flag of r.in.gdal, no ?

Moritz,

The problem is that I don't know the region's extents now. The project
location is a point and I have bathymetric and topographic maps that extend
well behond that point.

Thanks,

Rich