[GRASS-user] r.input and d.rast issues

rshepard · September 28, 2021, 9:13pm

Since r.in.gdal followed by r.proj are not working for me, despite success
with the rastertest location I applied Carlos' suggestion to use r.import.

Starting GRASS from /data/grassdata/project/topography/ I applied r.import
to the rastertest map, columbia_2010_e_dtm_35.tif with this result:

r.import in=~/projects/washington/nevins-dock/data/topography/columbia_2010_e_dtm_35.tif mem=6000 out=dtm35 resamp=lanczos_f ext=region resolution=region --o

Building topology for vector map
<tmp_r_import_region_salmo_7551@topography>...
Registering primitives...
Building areas...
  100%
Attaching islands...
  100%
Attaching centroids...
  100%
Selected PROJ pipeline:
+proj=pipeline +step +inv +proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5
+lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000 +y_0=0
+ellps=GRS80 +step +proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5
+lat_1=45.8333333333333 +lat_2=47.3333333333333 +x_0=500000 +y_0=0
+ellps=GRS80 +step +proj=unitconvert +xy_in=m +xy_out=us-ft
************************
Reprojecting primitives ...
Building topology for vector map
<tmp_r_import_region_salmo_7551@PERMANENT>...
Registering primitives...
Building areas...
  100%
Attaching islands...
  100%
Attaching centroids...
  100%
ERROR: Input raster does not overlap current computational region. Nothing
        to import.
ERROR: Unable to import GDAL dataset
        </home/rshepard/projects/washington/nevins-dock/data/topography/columbia_2010_e_dtm_35.tif>

I interpreted the manual page explanation of region to shift the output to
that region. Do I need to set the current region to the region of the map
I'm importing? How do I learn this before I import the .tif?

I set the region to the r.info map=local_depth
+-------------------- | Map: local_depth | Mapset: topography | Location: new_nevins_dock | DataBase: /data/grassdata | Title: | Timestamp: none |-------------------- | | Type of Map: raster | Data Type: FCELL | Rows: 15265 | Columns: 10768 | Total Cells: 164373520 | Projection: | N: | E: | Range of data: | | Data Description: | generated by r.proj | | Comments: | r.proj | pth" method="nearest" | +-------------------- raster map, local_depth:
--------------------------------------------------------+
Date: Fri Sep 24 15:49:01 2021 |
Login of Creator: rshepard |
|
|
|
|
--------------------------------------------------------|
|
Number of Categories: 0 |
Band reference: (none) |
|
|
|
NAD83(HARN) / Washington South |
106729 S: 60934 Res: 3 |
1544180 W: 1511876 Res: 3 |
min = NULL max = NULL |
|
|
|
|
|
-n location="project/ mapset="bathymetry" input="local_de\ |
memory=300 |
|
--------------------------------------------------------+

but when I try to display it now GRASS tells me it cannot find the map to
display:
d.rast map=local_depth
python3: can't open file '/data/grassdata/project/bathymetry/.tmp/salmo/MONITORS/wx0/render.py': [Errno 2] No such file or directory

What might be going on here?

TIA,

Rich

Micha_Silver1 · September 29, 2021, 8:50am

On 9/29/21 12:13 AM, Rich Shepard wrote:

Since r.in.gdal followed by r.proj are not working for me, despite success
with the rastertest location I applied Carlos' suggestion to use r.import.

Starting GRASS from /data/grassdata/project/topography/ I applied r.import
to the rastertest map, columbia_2010_e_dtm_35.tif with this result:

r.import in=~/projects/washington/nevins-dock/data/topography/columbia_2010_e_dtm_35.tif mem=6000 out=dtm35 resamp=lanczos_f ext=region resolution=region --o

Building topology for vector map
<tmp_r_import_region_salmo_7551@topography>...
Registering primitives...
Building areas...
100%
Attaching islands...
100%
Attaching centroids...
100%
Selected PROJ pipeline:
+proj=pipeline +step +inv +proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5
+lat_1=47.3333333333333 +lat_2=45.8333333333333 +x_0=500000 +y_0=0
+ellps=GRS80 +step +proj=lcc +lat_0=45.3333333333333 +lon_0=-120.5
+lat_1=45.8333333333333 +lat_2=47.3333333333333 +x_0=500000 +y_0=0
+ellps=GRS80 +step +proj=unitconvert +xy_in=m +xy_out=us-ft
************************
Reprojecting primitives ...
Building topology for vector map
<tmp_r_import_region_salmo_7551@PERMANENT>...
Registering primitives...
Building areas...
100%
Attaching islands...
100%
Attaching centroids...
100%
ERROR: Input raster does not overlap current computational region. Nothing
to import.
ERROR: Unable to import GDAL dataset
</home/rshepard/projects/washington/nevins-dock/data/topography/columbia_2010_e_dtm_35.tif>

I interpreted the manual page explanation of region to shift the output to
that region. Do I need to set the current region to the region of the map
I'm importing? How do I learn this before I import the .tif?

I set the region to the raster map, local_depth:
r.info map=local_depth
+----------------------------------------------------------------------------+

| Map: local_depth Date: Fri Sep 24 15:49:01 2021 |
| Mapset: topography Login of Creator: rshepard |
| Location: new_nevins_dock |
| DataBase: /data/grassdata |
| Title: |
| Timestamp: none |
|----------------------------------------------------------------------------|

| |
| Type of Map: raster Number of Categories: 0 |
| Data Type: FCELL Band reference: (none) |
| Rows: 15265 |
| Columns: 10768 |
| Total Cells: 164373520 |
| Projection: NAD83(HARN) / Washington South |
| N: 106729 S: 60934 Res: 3 |
| E: 1544180 W: 1511876 Res: 3 |
| Range of data: min = NULL max = NULL |
| |
| Data Description: |
| generated by r.proj |
| |
| Comments: |
| r.proj -n location="project/ mapset="bathymetry" input="local_de\ |
| pth" method="nearest" memory=300 |
| |
+----------------------------------------------------------------------------+

but when I try to display it now GRASS tells me it cannot find the map to
display:
d.rast map=local_depth
python3: can't open file '/data/grassdata/project/bathymetry/.tmp/salmo/MONITORS/wx0/render.py': [Errno 2] No such file or directory

What might be going on here?

I think there are two things going on here:

1. You reprojected the `local_depth` raster from the topography location to new_nevins_dock, but without setting the target mapset in advance. So the resulting local_depth had NULL values everywhere. Regarding projecting with r.proj, note this comment in the man page:

/"To avoid excessive time consumption when reprojecting a map the region and resolution of the target location should be set appropriately beforehand./

/A simple way to do this is to check the projected bounds of the input map in the current location's projection using the////*-p*////flag. The////*-g*////flag reports the same thing, but in a form which can be directly cut and pasted into a//g.region <https://grass.osgeo.org/grass78/manuals/g.region.html>command\. After setting the region in that way you might check the cell resolution with "g.region -p" then snap it to a regular grid withg.region <https://grass.osgeo.org/grass78/manuals/g.region.html>'s\*\-a\*flag\. E.g.g.region -a res=5 -p. Note that this is just a rough guide./

/A more involved, but more accurate, way to do this is to generate a vector "box" map of the region in the source location using//v.in.region -d <https://grass.osgeo.org/grass78/manuals/v.in.region.html>\. This "box" map is then reprojected into the target location withv.proj <https://grass.osgeo.org/grass78/manuals/v.proj.html>\. Next the region in the target location is set to the extent of the new vector map withg.region <https://grass.osgeo.org/grass78/manuals/g.region.html>along with the desired raster resolution (g.region -mcan be used in Latitude/Longitude locations to measure the geodetic length of a pixel).r.projis then run for the raster map the user wants to reproject. In this case a little preparation goes a long way."/

2. Then you used to `extent=` parameter to import columbia_2010_e_dtm_35.tif only within the current computational region, which you had set to match the raster `local_depth`. However it seems that these do not overlap at all, so the import failed.

What does:

gdalinfo ~/projects/washington/nevins-dock/data/topography/columbia_2010_e_dtm_35.tif | grep -A4 "Coordinates"

show??

Does that help?

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

rshepard · September 29, 2021, 12:51pm

On Wed, 29 Sep 2021, Micha Silver wrote:

Micha,

1. You reprojected the `local_depth` raster from the topography location
to new_nevins_dock, but without setting the target mapset in advance. So
the resulting local_depth had NULL values everywhere. Regarding projecting
with r.proj, note this comment in the man page:

The 'local_depth' map is in the bathymetry mapset and I set the target
region to it for all reprojections.

What does:

gdalinfo ~/projects/washington/nevins-dock/data/topography/columbia_2010_e_dtm_35.tif | grep -A4 "Coordinates"

show??

It shows corners in US-feet so I cannot compare them to the corners of the
region which are in meters.

I've not before had such issues importing and reprojecting raster data.
Rather than continuing to try to push water uphill I've decided on an
alternate way of presenting the bathymetric information in the reports.

Thanks again,

Rich