GRASS GSoC 2012 High level map interaction

(See also other GRASS GSoC 2012 projects)

Student Name:	Pietro Zambelli, University of Trento, Italy
Organization:	OSGeo - Open Source Geospatial Foundation
Mentor Name:	Mentor: Sören Gebbert, Backup mentors: Luca Delucchi, Martin Landa
Title:	Proposal for a high level maps interaction

Proposal for a high level maps interaction

The idea is: extend the python GRASS API to make it more pythonic :-)

As you can see in the subsequent tests, a lot of functionality are redundant, similarly with `numpy` library I think that could be a good idea binding all the grass functions and add some of them as methods of the class.

For example you can use both:

    # import the raster function
    >>> import grass.obj.function.raster as r

    # import the raster object library
    >>> import grass.obj.raster as raster

    # initialize a RasterObject present in the mapset with
    >>> dtm = raster.Raster('dtm10')

    # then you can print the info using the class method
    >>> print(dtm.info())

    # or using grass function
    >>> print(r.info(dtm))

Grass as python library

I would like to have the ability to call grass as a normal python library outside grass environment::

    # import the library
    >>> import grass

    # specify the mapset
    >>> grass.obj.config(mapset='/path/to/gisdata/location/mapset')

Region as an object

Region Attributes, read and write

    # import the new virtual extension of the GRASS API 
    >>> import grass.obj.region as region

    # save as default region option: `-s`
    >>> region.default = True

    # see region parameters as attribute of the class
    >>> region.nsres
    10

    # change the region parameters in a more pythonic way
    >>> region.nsres = 20

    # view the region extension
    >>> region.bbox
    [(4928030,589980),(4913690, 609000)]

    # change the region extension directly
    >>> region.bbox = [(4928030,589980),(4913690, 609000)]

    # view the convergence angle (degrees CCW)
    >>> region.conv_angle
    0.0

    >>> region.name = 'Trento'

Region Methods

    # print the current region `-p`
    >>> print(region)  # using the `__str__` method

    # set as default region `-s`
    >>> region.set_as_default() # It's equal to set `region.default = True`

    # set from default region `-d`
    >>> region.set_from(regionname='default')

    # set region from a map
    >>> region.set_from_map(MapObject)

    # align region to resolution (default = align to bounds, works only for 2D resolution)
    >>> region.align()

    # adjust region cells to cleanly align with this raster map
    >>> region.align(RastObj)

    # shrink region until it meets non-NULL data from this raster map
    >>> region.zoom(RastObj)

    # save current region settings in named region file
    >>> region.save('/your/path') # will produce a file: '/your/path/%s.region' % region.name

Maps as objects

MapObject class

Map Attributes

Similarly with the region the MapObject is characterize by some attributes, but only some of these attributes are modifiable.

Read and Write Attributes

Suppose that we initialize a MapObject called `dtm`

    # import the raster object library
    >>> import grass.obj.raster as raster
    
    # import a raster map from file into the Location and load
    >>> dtm = raster.import(fpath = '/your/data/DTM10.txt', name = 'dtm10')

    # view the Map name used in the mapset
    >>> dtm.name
    'dtm10'

    # to rename the map you have just to reassign the Map attribute
    >>> dtm.name = 'dtm__10x10m'

    >>> dtm.title
    "Digital Terrain Model compute from LIDAR2008"

    # to change the title, as before, you can just reassign the Map attribute
    >>> dtm.title = 'Digital Terrain Model compute from LIDAR2008; Licence: Creative Common 0"

    # to see in witch mapset the map is
    >>>,dtm.mapset
    '/your/gisbase/location0/mapset'

    # if you want to copy the map in another mapset changing the projection
    >>> dtm.mapset = '/your/gisbase/location1/mapset

Read only attributes

This attributes are only readable, and are relative to the attributes of the map and not to the region.

    >>> dtm.north
    4928030
    >>> dtm.south
    4913690
    >>> dtm.east
    609000
    >>> dtm.west
    589980
    
    # what append if the user try to change the resolution of the map
    >>> dtm.north = 4928000
    SyntaxError: invalid syntax.
        you can not modify the map attributes,
        please consider to change region parameters.

    >>> dtm.bbox
    [(4928030,589980),(4913690, 609000)]
    >>> dtm.type
    'raster'
    >>> dtm.creator
    'pietro'
    >>> dtm.creationtime
    datetime.datetime(2012, 3, 9, 17, 25, 6, 406987)

Map Methods

    # return the args and the kargs that generated the map
    >>> dtm.generated_by()
    (["r.in.ogr", "bla", "bla"], {'overwrite' : True})

    # return the bash command string that generate the map
    >>> dtm.bash()
    "r.in.ogr bla bla -overwrite=True"

    # to rename the map in the mapset
    >>> dtm.rename('newname0')

    # to copy the mapname
    >>> dtm.copy('newmap1')

    # to copy the mapname in an other mapset
    >>> dtm.copy('newmap2', mapset='path/to/your/mapset')

Raster Map

Raster attributes

Raster class is a child class of the MapObject and inherits all the attributes and methods, and add some more specific attributes and methods

Read only attributes

    >>> dtm.min
    0.000000
    >>> dtm.max
    52.520164
    >>> dtm.nsres
    30
    >>> dtm.ewres
    30
    >>> dtm.datatype
    FCELL
    >>> dtm.rows
    477
    >>> dtm.cols
    634
    >>> dtm.cells
    302418

Raster Methods

One thing that some time could be useful is the access to the low level of C library from python, This is a stupid example, and for this kind of things is better to use the Map algebra, but I think that sometimes have an access to a low level function is useful.

    # import the raster object library
    >>> import grass.obj.raster as raster

    # initialize an empty raster map, copy from another
    >>> newmap = raster.Raster('newmap')
    >>> for row in dtm:
    ...     for pixel_coords in row.not_null():
    ...         # pixel_coords is a couple of values, with the row and col
    ...         # do something...

Multilayer raster (i.group)

    # initialize a multi raster map
    >>> multi = raster.MultiRaster('multi')

    # add layers
    >>> multi['elev'] = dtm
    >>> multi['slope'] = slope_aspect(dtm, slope = 'slope')

    # query to one point
    >>> multi[10,10]
    {'elev': 100, 'slope': 30}
    >>> multi['elev'][10,10]
    100
    >>> multi.layers() or multi.keys()
    ('elev', 'slope')
    >>> multi[10,10].items()
    (val0, val1)

Time raster

I know that there are a lot of new staff regarding the time maps, but I don't know exactly how they works therefore any hint are welcome.

Vector MapObject

For the vector we can add some particular attributes and methods.

    # import the vector object library
    >>> import grass.obj.vector as vector
    
    # import a vector map from file into the location and load
    >>> roads = vector.import(fname = 'roads.shp', name = 'roads')

    # or if the map already exist in the mapset
    >>> roads = vector.Vector('roads')

    # add a category column
    >>> roads.add_column('length')
    >>> for road in roads:
    ...     road.cols.length = road.length()

    >>> for road in roads:
    ...     for point in road:
    ...         # do something with point

Execute function

    >>> import grass.obj.function.vector as v
    >>> import grass.obj.function.raster as r
    >>> from math import pi

    >>> road_buff_50 = v.buffer(roads, dist = 50.)

    >>> slope_deg, aspect = r.slope_aspect(dtm, slope = 'slope10_deg',
    ...                                    aspect = 'aspect10')
    >>> slope_rad = pi / 180 * slope_deg

    # we can set that default values name are slope = 'slope' and aspect = 'aspect'
    # PROBLEM: some function like r.slope.aspect, may return a different number of maps
    # the function generate these map only if the name is given, could be?
    # we need to generate a python function that read all the C header and generate
    # the code, when the user build the package, using the GCC-XML?
    # idea? hint?
    >>> slope, aspect, dxx = r.slope_aspect(dtm,  slope = 'slope10_deg',
    ...                                     aspect = 'aspect10',
    ...                                     dxx = 'dxx')

    # of course we need to handle wrong input map
    >>> slope, aspect = r.slope_aspect(road)
    TypeMapError: The function: "slope_aspect" require a raster map.

Algebra and logic function

Use mathematic functions and map algebra.

    >>> import numpy as np
    >>> drop = dtm.nsres * np.sin(slope_rad)
    >>> dist = dtm.nsres * np.cos(slope_rad)
    >>> length = np.sqrt(drop**2 + dist**2)
    
    # use logic ::
    >>> slope30 = slope if slope<=30 else None
    
    # and similarly mixing raster and vector maps ::
    >>> dtmroad = dtm if road else None
    >>> dtmroad100 = dtm if road.cols.length <= 100 else None

    # or you can make a mask, with boolean values: 
    >>> slope30 = slope <= 30

    # or query a vector map with:
    >>> track = road.cols.highway == 'track'

Save and Export

    # save all the maps created during the session ::
    >>> grass.obj.session.save()

    # or save or export, just one map: ::
    >>> slope30.save()

    # or save specify another location and mapset, with: ::
    >>> slope30.save(location='locationname', mapset='mapsetname')

    # export:
    >>> slope30.export('filename', format='asciigrid')

    # allow to convert a raster map to vector, like
    >>> slope30.export('filename', format='shp')
    You are converting a raster map to vector.
    Choose your Feature type [POINT/line/area]: point

    # or given directly the geometry type
    >>> slope30.export('filename', format='shp', feature_type = 'point')