Loading ADF file from National Elevation Dataset into GRASS with Python?

Question

I'm trying to load data from the National Elevation Dataset into GRASS/Python, using this Python script (outside of GRASS):

# Imports
import os
import grass_session
from grass.script import setup as gsetup
import grass.script as gscript
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import imagery as i
from grass.pygrass.modules.shortcuts import display as d

# Settings
gisdb = "/home/user/research/grassdata"
location = "USA"
mapset = "PERMANENT"

# Key directories
project = "/home/user/research/radio/analysis/gis"

# Initialization
gisbase = os.environ['GISBASE']
gsetup.init(gisbase, gisdb, location, mapset)


# Sample GIS data
adf_file = os.path.join(project, "sample", "grdn48w095_13", "w001001.adf")
output_file = os.path.join(project, "output", "sample")

print("Importing raster data... ", end = "")
r.in_gdal(input = adf_file, output = output_file, flags = "e", overwrite = False)
print("Done.")

But I get the following error:

Traceback (most recent call last):

  File "/home/user/apps/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py", line 2961, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)

  File "<ipython-input-7-77fc8fd256bb>", line 31, in <module>
    r.in_gdal(input = adf_file, output = output_file, flags = "e", overwrite = False)

  File "/usr/lib/grass72/etc/python/grass/pygrass/modules/shortcuts.py", line 47, in __getattr__
    return self.cls('%s.%s' % (self.prefix, name.replace('_', '.')))

  File "/usr/lib/grass72/etc/python/grass/pygrass/modules/interface/module.py", line 475, in __init__
    tree = fromstring(self.xml)

  File "/home/user/apps/anaconda3/lib/python3.6/xml/etree/ElementTree.py", line 1315, in XML
    return parser.close()

  File "<string>", line unknown
ParseError: no element found: line 1, column 0

As far as I can tell, I'm following the documentation correctly in my startup script, and my simple code roughly follows the GRASS wiki article on python.

What am I doing wrong? I've never used GRASS or GIS software before, so maybe I'm missing something obvious here. I'm using Grass 7.2 on a Debian Linux system with Python 3 (installed with Anaconda).

Answer 1

I think your problem is the output string. When you import a raster layer into GRASS you are creating a GRASS raster in the GRASS database, and with the internal format used by GRASS, not a file in some other directory. So the output parameter should be (in your case) "sample". Without os.path.join(). When you run the command, that will create a raster in the USA location, and in the PERMANENT mapset named sample.

I would suggest you start GRASS itself, not thru a python script and run this command in the interface (either GUI or command line) to learn how it works. Then, once you feel comfortable with GRASS, move to python scripting.