2

I am making a Python script for processing my data.

import grass.script as grass
import numpy as np

def main():
     f = open("D:\\Glaciers_selected.csv", 'r')
     fl = f.readlines()
     for x in range(1, len(fl)):
           nf = fl[x].split(';')[1]

           #### Problem in next executable line
           ## I need to create a new location, using the projection from the nf file. 
           grass.read_command('grass76' flags='c' nf ./UTM60S/     ## Here I am stuck

           ## The rest of the script works fine
           grass.read_command('r.import', input= 'D:\\GIS\Farinotti_glaciers\\' + nf + '.tif ', output='RGI60_11', overwrite=True)
           grass.read_command('g.region', raster='RGI60_11', flags='a', overwrite=True) 
           ...
           ...

if __name__ == "__main__":
    main()

My script runs fine if I start GRASS, add a new location with the location wizard and then run each line, one by one from the Python tab in the layer manager. However, I do not know how to add a location once I am inside a loop of a Python script.

When I run grass76 -c from the console, the location Wizard appears ...

I read Creating a GRASS location from a bash script and there seems to be a way of doing this, but I was not able to follow the explanation (which is for BASH, by the way, not for Python).

2

What you need is the g.proj command with the -c flag. However without knowing what's in your Glaciers_selected.csv file, I can't know how to specify the coordinate reference system. The -c flag for g.proj (just like grass -c ...) requires either a georeferenced data file, an EPSG code, or the full proj.4 CRS specification.

Here's an example. I start in a location defined as Pseudo Mercator, EPSG 3857, and create a new location based on a georeferenced shapefile in EPSG 4326. I switch to that location, then import the layer projected in that CRS (EPSG 4326):

micha@RMS World $ ipython
Python 2.7.15+ (default, Nov 27 2018, 23:36:35) 
Type "copyright", "credits" or "license" for more information.

IPython 5.5.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

import grass.script as gscript

# Show current location CRS
gscript.run_command('g.proj', flags="p")
-PROJ_INFO-------------------------------------------------
name       : WGS 84 / Pseudo-Mercator
a          : 6378137
es         : 0
proj       : merc
lat_ts     : 0.0
lon_0      : 0.0
x_0        : 0.0
y_0        : 0
k          : 1.0
wktext     : defined
no_defs    : defined
-PROJ_EPSG-------------------------------------------------
epsg       : 3857
-PROJ_UNITS------------------------------------------------
unit       : Meter
units      : Meters
meters     : 1

# Now create a new location, using a georeferenced file
gscript.run_command('g.proj', flags="c", georef="Mediterranean.shp", location="longlat")
Location <longlat> created
You can switch to the new location by
`g.mapset mapset=PERMANENT location=longlat`

# Change to that new location, and display new CRS parameters
gscript.run_command('g.mapset', location="longlat", mapset="PERMANENT")
Mapset switched. Your shell continues to use the history for the old mapset
You can switch the history by commands:
history -w; history -r
/home/micha/work/tmp/longlat/PERMANENT/.bash_history;
HISTFILE=/home/micha/work/tmp/longlat/PERMANENT/.bash_history

gscript.run_command('g.proj', flags="p")
-PROJ_INFO-------------------------------------------------
name       : WGS 84
datum      : wgs84
ellps      : wgs84
proj       : ll
no_defs    : defined
-PROJ_UNITS------------------------------------------------
unit       : degree
units      : degrees
meters     : 1.0

gscript.run_command('v.import', input_="Mediterranean.shp", output="med")
Check if OGR layer <Mediterranean> contains polygons...
 100%
Creating attribute table for layer <Mediterranean>...
Default driver / database set to:
driver: sqlite
database: $GISDBASE/$LOCATION_NAME/$MAPSET/sqlite/sqlite.db
Importing 2 features (OGR layer <Mediterranean>)...
 100%
# ....additional output lines of cleaning etc...
Input <Mediterranean.shp> successfully imported without reprojection

Note that the new location has the CRS of the input file, and so the v.import command does NOT need the -o parameter.

You can, of course, wrap the above in a loop, and create a new location for each input file if needed. But you must have the correct CRS available either as a georeferenced file, EPSG code, or the full specification.

0

This is how you create a new location from inside a python script:

grass.core.create_location(dbase='C:\\grassdata\\Glaciers\\PERMANENT\\sqlite', location='Lnew', filename='D:\\Farinotti_glaciers\\RGI60-11.00015.tif', overwrite=True)

dbase: is the path to GRASS database
location: location name to create
file: if given create new location based on georeferenced file

Additional arguments can be found in: Python documentation - script package - submodules

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.