Python scripting in GRASS

Question

A script that worked fine one week ago, has errors and warnings now. I must have changed something.

import grass.script as grass
import numpy as np

def main():
grass.read_command('r.import',input='D:\\dhm75_lcc.img',output='dhm75_lcc', overwrite=True)          
f = open("D:\\Glaciers_RGI_60_selected_classified.csv", 'r')
fl = f.readlines()
for x in range(1, len(fl)):
     nf = fl[x].split(';')[1]     # Get name of file to import
     grass.run_command('g.proj', flags='c', georef='D:\\Farinotti_composites\\' + nf + '.tif', location="test")
     grass.run_command('g.mapset', location="test", mapset="PERMANENT", overwrite=True)
     grass.read_command('r.import', input='D:\\Farinotti_composites\\' + nf + '.tif',output='RGI60_11',overwrite=True)
     grass.read_command('r.reclass',input='dhm75_lcc',output='dhm75_lcc_rc',rules='D:\\Elev_Reclass_rules.txt',overwrite=True)
     grass.read_command('r.reclass',input='RGI60_11',output= 'RGI60_11_rc',rules='D:\\Thickness_Reclass_rules.txt',overwrite=True)
     rs = grass.read_command('r.stats', flags='c', input = ['RGI60_11_rc', 'dhm75_lcc_rc'], separator = ' ')
     np.save('D:\\Cell_count_composite_'+nf, rs)

if __name__ == "__main__":
main()

This is the output I get:

WARNING: Projection of dataset does not appear to match current location.

Location PROJ_INFO is:
name: WGS 84 / UTM zone 33N
datum: wgs84
ellps: wgs84
proj: utm
zone: 33
no_defs: defined

Dataset PROJ_INFO is:
name: Lambert Conformal Conic
datum: wgs84
ellps: wgs84
proj: lcc
lat_1: 46
lat_2: 49
lat_0: 47.5
lon_0: 13.3333333333
x_0: 400000
y_0: 400000
towgs84: 0,0,0,-0,-0,-0,0
no_defs: defined

ERROR: proj
Importing raster map <dhm75_lcc>...
Estimated target resolution for input band <dhm75_lcc>: 73.3401507889
Reprojecting <dhm75_lcc>...
ERROR: Unable to create location <test>: File exists

As the error seems to be that I cannot overwrite the location, I tried to make a new location for each file:
So I changed two lines of the script: instead of using location="test", I used location="test"+str(x). As x differs for each file I get different location names for each file, so there is no need for overwriting.

...
for x in range(1, len(fl)):
  nf = fl[x].split(';')[1]           # Get name of file to import
  grass.run_command('g.proj', flags='c', georef='D:\\Farinotti_composites\\' + nf + '.tif', location="test"+str(x))
  grass.run_command('g.mapset', location="test"+str(x), mapset="PERMANENT", overwrite=True)
...

Now I get this error:

Importing raster map <dhm75_lcc>...
Estimated target resolution for input band <dhm75_lcc>: 73.3401507889
Reprojecting <dhm75_lcc>...
Location <test1> created
You can switch to the new location by
`g.mapset mapset=PERMANENT location=tets1`
WARNING: Concurrent mapset locking is not supported on Windows
Mapset switched. Your shell continues to use the history for the old mapset
Importing raster map <RGI60_11>... 
ERROR: Raster map <dhm75_lcc> not found

What I did now was to move the line where I read the DEM (dhm75_lcc) below the line where I switch to the new location. (Note: as this file is large and takes some time to import, I would prefer to read it just once at the beginning of the script and not in each loop):

...
def main():
    f = open("D:\\Glaciers_RGI_60_selected_classified.csv", 'r')
    fl = f.readlines()
    for x in range(1, len(fl)):
        nf = fl[x].split(';')[1]     # Get name of file to import
        grass.run_command('g.proj', flags='c', georef='D:\\Farinotti_composites\\' + nf + '.tif', location="test"+str(x))
        grass.run_command('g.mapset', location="test"+str(x), mapset="PERMANENT", overwrite=True)
        grass.read_command('r.import',input='D:\\dhm75_lcc.img',output='dhm75_lcc')  
...                 

The script runs now and produces a result, but (i) there are two warnings, (ii) one strange message that I was not getting before: "0.600000 rounded up to 0" and (iii) the results differ (although not much) from the results I got last week, when there were no warnings.

ERROR: proj
Importing raster map <dhm75_lcc>...
Estimated target resolution for input band <dhm75_lcc>: 73.3401507889
Reprojecting <dhm75_lcc>...
Importing raster map <RGI60_11>...
0.600000 rounded up to 0
Location <tett3> created
You can switch to the new location by
`g.mapset mapset=PERMANENT location=tett3`
WARNING: Concurrent mapset locking is not supported on Windows
Mapset switched. Your shell continues to use the history for the old mapset
WARNING: Projection of dataset does not appear to match current location.

So, why does the projection of dataset does not appear to match current location? Do i need to do something about "Concurrent mapset locking is not supported on Windows"?

Answer 1

In order to do any processing on rasters from different CRS's you'll have to reproject to a common CRS. The common CRS in your case might be UTM 33 since that's the projection of the DEM. However does your DEM, in UTM 33, cover the area of UTM31 - UTM 34 ??. I'm not sure that reprojecting from UTM 31 to UTM 33 would be appropriate - you'll get some visible distortion. You might choose some other projection that can cover that span from west to east.

Having said that, in general what you have to do is read thru the CSV line by line, as you have done, and for each line:

1. create a new temporary location for each input raster. This process is explained quite well in: this GRASS wiki page.
2. Change to that location, and import the raster. (Save the temp location name in a python variable, usinggrass.gisenv['LOCATION'] and grass.gisenv['MAPSET'] )
3. Change back to your common projection, and project the new raster from the temp location into your common location. See r.proj for details.
4. Now destroy the temporary location, and go to the next input line

WHen the loop finishes you can do any further processing you need on the rasters that should all now be in one uniform CRS under one location/mapset