Creating new GRASS locations while looping in Python script

Question

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).

Answer 1

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.

Answer 2

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