0

I want to crop a raster based on another one, inside a loop. I created a mask (MASK) and then cropped the old raster (old), creating a new raster (new).

r.mask vector=EU@PERMANENT #my mask 
r.mapcalc "new = old" 

For a single raster, it works. However, I have many old rasters, and I want to create a loop. I tried this:

for i in `cat /myPath/myRasters.txt | tr -d '\r'`
do
echo SET REGION 
g.region vector=EU@PERMANENT
echo CREATE A MASK 
r.mask --o vector=EU@PERMANENT  
echo CROP $i RASTER 
r.mapcalc "$i'_new' = $i'_old'"
r.mask -r 
done 

I can see the problem is in the syntax used, as I get this ERROR message after echo CROP $i RASTER:

syntax error, unexpected STRING, expecting '('

I tried also to use r.mapcalc with:

"$i'_new' = $i'_old'"
'$i'_new' = $i'_old''
"${i}'_new' = ${i}'_old'"
"$i'_new' = if(MASK,$i'_old')"

this last one following r.mapcalc manual page ( https://grass.osgeo.org/grass82/manuals/r.mapcalc.html ), where it says

To mask raster map layer b by raster map layer a:
result = if(a,b)

and this: https://grasswiki.osgeo.org/wiki/Cropping_maps

But nothing works. Previous similar questions (e.g., Loop using variable - r.mapcalc GRASS GIS ) didnt' help. I am using GRASS 7.8.6 on Ubuntu 18.04.6 (from remote, I have Windows 10).

6
  • First, do you really have a GRASS vector named "EU@PERMANENT" ?? That's a strange name for a vector! In GRASS the word PERMANENT is reserved for the default mapset in a location. I would avoid naming GRASS maps like that.
    – Micha
    Commented Jul 9, 2022 at 18:25
  • Second, I think this can be easier than you have made it. Can we see the contents of /myPath/myRasters.txt??
    – Micha
    Commented Jul 9, 2022 at 18:26
  • Dear Micha, no, the vector name's just EU, I can also write vector=EU, correct? The file myRasters contains, in each row, the name of each raster as it is loaded in my current mapset (e.g., raster_1, raster_2, and so on), while the file itself would be raster_1.tiff (or raster_1.img), and so on. The file myRasters.txt has correct endings (I checked it with cat -v) and has a last blank row.
    – LT17
    Commented Jul 10, 2022 at 12:06
  • Be sure to quote the echo commands: echo CROP $i RASTER --> echo "CROP $i RASTER".
    – markusN
    Commented Jul 10, 2022 at 22:18
  • Also: r.mapcalc "$i'_new' = $i'_old'" --> r.mapcalc "${i}_new = ${i}_old" (the curly braces are enough, no need for single quotes).
    – markusN
    Commented Jul 10, 2022 at 22:20

1 Answer 1

0

Following your answer in the comments, you have a text file with the raster names listed, row by row. So:

# Set region and mask once. No need to do this in the loop
echo "SET REGION" 
g.region vector=EU
echo "CREATE A MASK"
r.mask --o vector=EU

while read i; echo "CROP ${i}"; r.mapcalc "${i}_new = ${i}"; done < /myPath/myRasters.txt

r.mask -r   

That should do it.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.