1

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"?

  • In order to help, I think we need some more details: I understand that the CSV file contains a list of paths to Geotiff rasters? Are each of those in a different projections? Are they projected differently from the DEM? If you explain what you want to accomplish, rather than how you are trying to do it, it might be easier to make suggestions. – Micha Jul 21 at 17:43
  • Are you trying to import rasters that are in different projection than your created location? I would try to import rasters in one location rather than creating separate locations for each of them. Have you checked if the raster is a valid one? If any doubt on the raster import, I try to check by manually importing using 'r.import'. I ditto on the point that we need a bit more help on what you are trying to achieve. – Hasnein Tareque Jul 22 at 7:09
  • Yes, the CSV file has the path to each raster I want to process. The projection of the DEM is UTM 33, while the rasters I want to process have the projection UTM 31, UTM 32, UTM 33 or UTM 34. I am happy to provide more details if needed, I just do not really know what could be relevant. – Iced_ Jul 25 at 3:04
  • I have tested the rasters. They are OK, and I can process the files one by one using the GRASS GUI. I just have problems with the script. – Iced_ Jul 25 at 3:09
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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