0

I have a remote machine where my GRASS binaries (7.8) are available and I want to apply this solution https://gis.stackexchange.com/a/255899/36712 to my data set of thousands of small images. In most PyGRASS tutorials I find the API works more like an remote to the GUI, in other words you apply functions to layers present in the GrassGUI for example gscript.parse_command('r.univar', map='elevation', flags='g') where I apply the function r.univar to the layer elevation I allready found the grass-session libary allowing me to create a session without "starting" GRASS. Most of the tutorials I find are also either old referring to Python 2.7 or don't realy cover how to load a dataset or what would even better how I could pass a variable containing the file object to GrassGIS. How I imaging the code should look like is this:

from grass_session import Session

with Session(gisdb="/tmp", location="location", create_opts="EPSG:4326"):
    for png in list_of_png_files:
        geo_tiff = transform_png_to_geotiff_with_gdal(png)

        #either i need to do this:
        Loaded_geotiff = load_geotiff_with_grass(geo_tiff)

        #or I can work directly with the geo_tiff Variable
        output = processing.runalg('grass:r.mapcalculator',geo_tiff)
        thined = processing.runalg('grass7:r.thin',output)
        lines = processing.runalg('grass7:r.to.vect',thined)

        save_lines_to_shp(lines)

How can I achieve this? The PyGRASS API seems so counterintuitive and it's hard to understand how this is supposed to work.

  • I know tried with gs.run_command('r.external',input="./data/gtifs/13.tif",output="layer_1") but i just get an error CalledModuleError: Module run None r.external input=./data/gtifs/13.tif output=layer_1 ended with error Process ended with non-zero return code 1. See errors in the (error) output. – Thagor Apr 22 at 14:39
0

here is the answer for what i wanted to do

import tempfile
import glob
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import vector as v
from grass_session import Session

# glob.iglob gives me an glob itterator without reading the whole dir first
for i,png in enumerate(glob.iglob('/path_to_tiles/*/*.png'), recursive=False)):

    # creates the geotiff as an temp file to reduce the data waste
    with tempfile.NamedTemporaryFile(dir="./data/tmp_dir/") as tmp_geo_tif:

        # calls the function that transforms the png to a geotiff
        transform_to_geotiff(png,tmp_geo_tif.name)

        # because no grassgis instnace is runnig we must create one this is done with the sessions libary
        with Session(gisdb="./tmp", location="location", create_opts="EPSG:4326"):

            # we allways need a internal representation of the data_set we are working with 
            # here in the first step it is layer_1
            curr_raster='layer_1'

            # r. links to the shortcuts and this one links to the external raster import libary
            # it is opening the temp file containing the geotiff
            r.external(input=tmp_geo_tif.name, output=curr_raster, overwrite=True)

            # for mapcalc to work properly we need to set the region to the region of the active layer
            g.region(raster=curr_raster)

            # this results in the layer 'mapcalc' containing either the value of curr_raster or no data 
            # if the cell in curr_raster is bigger than 0 (transforms zero to nodata basicaly)
            r.mapcalc(f'mapcalc = if({curr_raster} > 0, layer_1, null())', overwrite = True)

            # r.thin fails if the raster after map calc only has no data cells thats why we need a try/except
            try:

                # r.thin takes the mapcalc and thins the raster lines to 1 pixel width
                r.thin(input='mapcalc', output='thinned', overwrite=True)

                # r.to_vect then transformst them into polylines
                r.to_vect(input='thinned', output='lines',type='line', overwrite = True)

                # lasty we save it to a gdpk
                output = f'./data/lines_gdpk/lines_{i}.gdpk'
                v.out_ogr(input='lines', output=output, overwrite = True)
            except:
                pass
| improve this answer | |

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.