Geomorphons for a geomoron: running grass command in external python

Question

I'm a fairly competent Linux and Python user and looking to utilize the GrassGIS geomorphon function (https://grass.osgeo.org/grass78/manuals/r.geomorphon.html) to parse out subaqueous features. I've never used GRASS and don't quite understand the architecture and work flow, so please bear with me.

In an ideal world I install the GRASS backend only via pip eg pip install pyGrass and call import pyGrass and pyGrass.r.geomorphon(arg1, arg2) (or something similar) similar to an interaction like is seen with the C based opencv library (https://opencv.org/) eg. pip install opencv then import cv2. My understanding is that not possible with the GRASS-GIS as it stands.

For example my envisioned workflow:

import pickle
from matplotlib import pyplot as plt 
url= "https://chldata.erdc.dren.mil/thredds/dodsC/frf/geomorphology/DEMs/surveyDEM/FRF_20191114_1177_FRF_NAVD88_LARC_GPS_UTC_v20191115_grid_latlon.nc"
ncfile = nc.Dataset(url)
lon = ncfile['longitude'][:]
lat = ncfile['latitude'][:]
elev = ncfile['elevation'][:].squeeze()

returnedArguments = pyGrass.r.geomorphon(elevation, arg1, arg2)
plt.figure()
plt.pcolormesh(lon, lat, elevation)
for idxType in returnedArguments:
    plt.pcolormesh(lon, lat, featureType[idxType])

I've found a few resources that I've begun cobbling together to work towards my end, none of which seem to illustrate my goals directly. It looks like most everything in the documentation is Python2 and I'm using Python3, CentOS

import os, sys, subprocess, pickle
url= "https://chldata.erdc.dren.mil/thredds/dodsC/frf/geomorphology/DEMs/surveyDEM/FRF_20191206_1179_FRF_NAVD88_LARC_GPS_UTC_v20191209_grid_latlon.nc"
ncfile = nc.Dataset(url)
lon = ncfile['longitude'][:]
lat = ncfile['latitude'][:]
elev = ncfile['elevation'][:].squeeze()

location_path = os.path.join(gisdb, location)
if os.path.exists(location_path):
    shutil.rmtree(location_path)
# Create new location from EPSG code, and exit
startcmd = ("%s -c %s -e %s" % (grass7bin, epsg_code, location_path))
print(startcmd)
p = subprocess.Popen(startcmd, shell=True,
                     stdout=subprocess.PIPE, stderr=subprocess.PIPE)
###########################
grassPath = "/home/spike/grassStuff"
#################################################
# from https://github.com/wenzeslaus/try-grass-in-jupyter/blob/master/notebook.ipynb
# create GRASS GIS session
gisbase = str(subprocess.check_output(["grass78", "--config", "path"]), 'utf-8').strip().split()[-1]
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))
# import GRASS GIS packages we need
import grass.script as gs
import grass.script.setup as gsetup
import grass.script.array as garray
rcfile = gsetup.init(gisbase, grassPath, location='demo')
# we want functions to raise exceptions and see standard output of the modules
gs.set_raise_on_error(True)
gs.set_capture_stderr(True)
# simply overwrite existing maps like we overwrite Python variable values
os.environ['GRASS_OVERWRITE'] = '1'
# from https://grasswiki.osgeo.org/wiki/GRASS_Python_Scripting_Library#Interfacing_with_NumPy
a = garray.array()   # I error out here, but this is not the focus of my question, i'm curious why
a.read(elevation)

As shown above, I would love to get here and is the main focus of my question:

# I'm guessing that I have to write the data to a raster
gs.grass.run_command(gscore.g.region raster=elevation -p)
# then gs.execute or some equivalent that returns data
returnedArgs = grass.run_command(r.geomorphon elevation=elevation forms='all') 
# then return values back as numpy array
dataIAmInterestedIn = a.write(returnedArgs)

but I'm quite sure that's totally wrong. What should I do instead?

Answer 1

I believe you're missing a few important things in your code. First, you need to define the GRASS region setting before trying to create a numpy array with garray.array(). Next the output raster from the r.geomorphons module is the string given to the forms parameter. You can then export that back to a numpy array (if needed).

Here's a full working example. I made this example by starting in a GRASS session. I clipped a small section of an elevation raster (from SRTM) then I made two additional rasters with the longitude and latitude values at each pixel - to mimic your situation. Then I started python from the GRASS session, and with the grass.script.array module, I made numpy arrays of each of those rasters, and saved as a dict to a pickle file.

After that preparation, the following python script reads in the python dict from the pickle file, and uses the numpy array to define region settings and create the GRASS elevation raster. The GRASS session is initialized within the script (as you did, more or less). Then the r.geomorphon module is called using grass.script.run_command(...), and the result is saved back to a pickle file.

import os
import sys
import subprocess
import pickle as pkl
import numpy as np

#######################################
# Get numpy arrays
#######################################
with open('data_list.pkl', 'rb') as fl:
    data = pkl.load(fl)
lon = data['lon']
lat = data['lat']
elev = data['elev']

#######################################
# Initialize GRASS session in temp directory
#######################################
# See: 
# https://grasswiki.osgeo.org/wiki/Working_with_GRASS_without_starting_it_explicitly#Python:_GRASS_GIS_7_without_existing_location_using_metadata_only
grass7bin = '/usr/bin/grass'
# Assume longitude/latitude location in WGS84 datum
epsg_code = 'EPSG:4326'
gisbase = str(subprocess.check_output([grass7bin, "--config", "path"]), 'utf-8').strip().split()[-1]
os.environ['GISBASE'] = gisbase
sys.path.append(os.path.join(gisbase, "etc", "python"))
# put temp GRASS location in current working directory
gisdb = './grassdata'
location = 'tmp_location'
mapset   = 'PERMANENT'
location_path = os.path.join(gisdb, location)
# Create new location from EPSG code, and exit 
startcmd = ("%s -c %s -e %s" % (grass7bin, epsg_code, location_path))
print(startcmd)
p = subprocess.Popen(startcmd, shell=True,
                     stdout=subprocess.PIPE, stderr=subprocess.PIPE)
out, err = p.communicate()
if p.returncode != 0:
    print('ERROR: %s' % err)
    print('ERROR: Cannot generate location (%s)' % startcmd)
    sys.exit(-1)
else:
    print('Created location %s' % location_path)
# Now the location with PERMANENT mapset exists.

#######################################
# Launch session and set region from the numpy array extents
#######################################
import grass.script as gs
from grass.script import array as garray
from grass.script import setup as gsetup
gsetup.init(gisbase, gisdb, location, mapset)

# Region settings
ncols = np.shape(elev)[1]
nrows = np.shape(elev)[0]
north = lat.max()
south = lat.min()
east = lon.max()
west = lon.min()
gs.run_command('g.region',
    n=north, s=south, e=east, w=west,
    rows=nrows, cols=ncols, flags="p")

# Create new array and assign values from the elev numpy array
elev_arr = garray.array()
elev_arr[:] = elev
# Save it to a GRASS raster.
# Careful not be be confused:
# both the original numpy array and the GRASS raster are named "elev"
elev_arr.write('elev', overwrite=True)
# Now run GRASS modules on that raster
gs.run_command('r.info', map_="elev")

#######################################
# Geomorphon
#######################################
# string variable to hold name of raster result
geomorph_result = 'geomorph_result'
gs.run_command('r.geomorphon', elevation="elev", forms=geomorph_result)
gs.run_command('r.category', map=geomorph_result)
# Convert GRASS raster to numpy array
geomorph_arr = garray.array()
geomorph_arr.read(geomorph_result)
geomorph_arr = np.asarray(geomorph_arr, dtype = np.integer)
with open("geomorphon_result.pkl", "wb") as fl:
    pkl.dump(geomorph_arr, fl)

I hope this helps you to move ahead.

I'm tempted to ask where the original numpy elevation array comes from. Why is this round trip from numpy to GRASS back to numpy necessary? Wouldn't it have been possible to import the elevation directly into GRASS??