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

Suppose that we initialize a MapObject called `dtm`

    # import the raster object library
    >>> import grass.obj.Raster as Raster
    
    # load an existing raster map, read mode
    >>> dtm = Raster(name = "dtm", mapset = "", mode = 'r', method = 'row')
    # create a new map, write mode
    >>> dtm2 = Raster(name = "dtm2", mapset = "", mode = 'w', method = 'row')
    # open the maps
    >>> dtm.open(); dtm2.open()

    # access row by row
    >>> for i in range(dtm.rows):
    ...     dtm2[i] = 2 * dtm[i] # dtm[i] return a numpy array

    # access cell by cell
    >>> for row in dtm:
    ...     for cell in row:
    ...         # do something

    # close
    >>> dtm.close(); dtm2.close()

    # open in read & write mode
    >>> dtm.open(mode = 'rw', method = 'segmentation')
    >>> for row in dtm:
    ...     row *= 2
    >>> dtm.close()

    # 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')

Weekly reports

Report0 - 2012-05-25

What did I do this week?

• studying the C-API of grass and ctypes;
• at the moment I'm in Prague at the Code Sprint with a lot of main developers and power users of the amazing GRASS community! and I have the opportunity to exchange ideas with them.
• start to develop a first implementation of the Raster and Region class.
• open a new website on | google code
• make a new repository, have a look: git clone https://code.google.com/p/pygrass/ or | download it.

What will I be working on next week?

Go on with the develop to open the Raster class using the read and write mode with the `row` method. Clean the Region class to use only ctypes and avoid the use grass.run_command. Continue to study the segment library of the C API.

Did I meet with any stumbling blocks?

I had some problem trying to use the segment library, but I hope to solve during the next week.

Report1 - 2012-06-01

What did I do this week?

• add more functionality to the raster object (rename, remove, del()):

# open an existing map
>>> elev = obj.Raster('elevation')
>>> elev.open()
# get the name of the map
>>> elev.name
'elevation'
# rename the map
>>> elev.name = 'elev'
>>> elev.name
'elev'
# or rename using the method
>>> elev.rename('elevation')
>>> elev.name
'elevation'
# we can check if the map exist with
>>> elev.exist()
True
# we can remove the map with
>>> elev.remove()
>>> elev.exist()
False
# or we can remove the map and de-instantiate the map object with
>>> del(elev)

• re-factoring some part of the code, the write row method has been changed from:

# open a new map
>>> new = obj.Raster('new', mode = 'w')
>>> new.open()

>>> c = 0
>>> for row in elev:
...     new[c] = row

when we are writing with 'row' method, we can not choose which row write, the above syntax could make the users confused, we just add one more row to the file, so I have changed the behavior adding a new method 'writerow', that is more explicit, so now we write:

>>> for row in elev: 
...     new.writerow(row)

• Try to add a new Row object, to extend the capabilities, at the moment the Raster object return a LP_c_float, that support slice, but not support iteration and basic algebra:

>>> type(elev[0])
<class 'grass.lib.ctypes_preamble.LP_c_float'>

# making a Row object
>>> row = Row(elev[0], elev._cols)
# add the support for a sum
>>> row2 = row + 2
# we can convert the row into an array
>>> row.to_array()
array([  1.03014536e+015,   1.06602391e+015,   1.12519191e+015, ...,
         0.00000000e+000,   0.00000000e+000,   4.94065646e-324])

• Continuing to study the segmentation library, add some new function and C struct, that should call by the Raster class using ctypes.
• Start to think how the map algebra should be implemented, may be I could use jinja to make a template to generate a C code that will be compile when user make GRASS, then the Raster class just call this function with ctypes.

What will I be working on next week?

Implement the read and write mode using the segmentation library.

Did I meet with any stumbling blocks?

I do not understand how should I extend the class object "LP_c_float", without write a new class (Row).

Report2 - 2012-06-08

What did I do this week?

• As suggested by the mentor I split the Raster class into:

   - RasterAbstractBase;
   - RasterRow;
   - RasterRowIO;
   - RasterSegment;
   - RasterNumpy;

• Discussed with my mentor the next steps;
• Studied the ctypes interface;
• Finished the module use to improve my C knowledge and to look the C API of grass | r.distdrop;
• Took a course at the university;

What will I be working on next week?

Develop the RasterSegment class.

Did I meet with any stumbling blocks?

Still problem with "LP_c_float"

# import and open a map
>>> import obj
>>> elev = obj.RasterRow('elevation')
>>> elev.open()

# instantiate a Row object
>>> row0 = elev[0]

# now if I look at the object return by Rast_get_row function
>>> elev._pbuf[:3]
[141.9961395263672, 141.2784881591797, 141.37904357910156]

# rather than the row is
>>> row0[:3]
Row([  1.03014536e+15,   1.06602391e+15,   1.12519191e+15])

as you see there is not correspondence, so I'm not convert correctly the LP_c_float into a buffer.

I create the row with:

def __getitem__(self, row):
    libraster.Rast_get_row(self._fd, self._pbuf, row, self._type)
    buf = np.core.multiarray.int_asbuffer( c.addressof(self._pbuf.contents),
                                           8*self._cols)
    self.row = Row( (self._cols, ), buffer = buf)
    return self.row

The good thing is that the elev._pbuf and row0 are looking at the same memory, if I change the first row value:

>>>> row0[0]
1030145362214880.4

>>> row0[0] = 10
>>> elev._pbuf[:3]
[0.0, 2.5625, 141.37904357910156]

>>> row0[0] = 1030145362214880.4
>>> elev._pbuf[:3]
[141.9961395263672, 141.2784881591797, 141.37904357910156]

if I change the second value:

>>> row0[1]                     
1066023908387873.1

>>> elev._pbuf[:4]
[141.9961395263672, 141.2784881591797, 141.37904357910156, 142.2982177734375]

>>> row0[1] = 10
>>> elev._pbuf[:4]
[141.9961395263672, 141.2784881591797, 0.0, 2.5625]

>>> row0[1] = 1066023908387873.1
>>> elev._pbuf[:4]
[141.9961395263672, 141.2784881591797, 141.37904357910156, 142.2982177734375

I don't understand this behavior. Any hint?

Report3 - 2012-06-15

What did I do this week?

Complete the Row class, the Row is an object that inherit from a | numpy.ndarray. It is possible to read and write using the RastrRow

>>> import obj
>>> elev = obj.RasterRow("elevation"); 
>>> elev.open()
>>> el_le_100 = obj.RasterRow('elev_le_100', 'w', mtype = 'CELL', overwrite = True)
>>> el_le_100.open() 
>>> for row in elev: 
...     el_le_100.put_row( row <= 100 )
... 

>>> el_le_100.close()

Finish the RasterSegment class:

>>> import obj
>>> elev = obj.RasterSegment('elevation')
>>> elev.open()
>>> for row in elev[:5]: print(row[:3])
[ 141.99613953  141.27848816  141.37904358]
[ 142.90461731  142.39450073  142.68611145]
[ 143.81854248  143.54707336  143.83972168]
[ 144.56524658  144.58493042  144.86477661]
[ 144.99488831  145.22894287  145.57142639]
>>> new = obj.RasterSegment('new')
>>> new.open()
>>> for irow in xrange(elev.rows): new[irow] = elev[irow] < 144
>>> for row in new[:5]: print(row[:3])
[1 1 1]
[1 1 1]
[1 1 1]
[0 0 0]
[0 0 0]
>>> elev[0, 0] == elev[0][0]
True
>>> new[0, 0] = 10
>>> new[0, 0]
10
>>> new.close()
>>> new.exist()
True
>>> new.remove()
>>> elev.close()

What will I be working on next week?

Next week I will start to work to integrate RowIO and Numpy. Start to add more documentation and tests and benchmarks.