GRASS GSoC 2012 High level map interaction

From GRASS-Wiki
Jump to navigation Jump to search

(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[1] 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')

Project plan

Period Task Status / Notes
May 25 Region class using g.region √
June 1 Implement Raster class √
June 8 Split Raster class, for each method √
June 15 Implement Buffer and RasterSegment classes √
June 22 Add documentation, doctest, and benchmark, RasterRowIO missing
June 29 Add RasterNumpy class missing
July 6 Add categories to Rasters classes missing
July 13 Implement module wrapper class or function generator missing
July 13 Mid-term evaluations deadline missing
July 20 Continue the module wrapper part missing
July 27 Implement Vector class missing
August 3 Implement Vector class missing
August 10 Implement Vector class missing
August 13 Suggested 'pencils down' date. Take a week to scrub code, write tests, improve documentation, etc. missing
August 17 Implement Vector class missing
August 20 Final evaluation

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.

Did I meet with any stumbling blocks?

No.

Report4 - 2012-06-22

What did I do this week?

Change names, and add new directories (obj => pygrass, docs, tests). Start the documentation of the | pygrass . Complete the RowIO class, and start to develop the RasterNumpy class. Start to develop tests using unittest. Benchmark the library, and compare the result with r.mapcalc. Met the co-mentors (Luca de Lucchi and Martin Landa) to illustrate the work done and discuss with them the next steps.

What can you do with the RowIO library: 

>>> import pygrass
>>> elev = pygrass.RasterRowIO('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]
>>> elev.close()


Start the RasterNumpy class: 

>>> import pygrass
>>> elev = pygrass.RasterNumpy('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]
>>> # in this case RasterNumpy is an extention of the numpy class
>>> # therefore you may use all the fancy things of numpy.
>>> elev[:5, :3]
RasterNumpy([[ 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]], dtype=float32)
>>> el = elev < 144
>>> el[:5, :3]
RasterNumpy([[ True,  True,  True],
   [ True,  True,  True],
   [ True,  True,  True],
   [False, False, False],
   [False, False, False]], dtype=bool)
>>> el.name = 'new'
>>> el.close()

# but I have problem to read the new map
>>> new = pygrass.RasterNumpy('new', mtype = 'CELL')
>>> new.open()
>>> new[:5,:3]
RasterNumpy([[1124990723, 1124943691, 1124950281],
       [1125050261, 1125016830, 1125035941],
       [1125110156, 1125092365, 1125111544],
       [1125159092, 1125160382, 1125178722],
       [1125187249, 1125202588, 1125225033]], dtype=int32)

What will I be working on next week?

Finish the RasterNumpy class, start to work on raster categories and may be history.

Did I meet with any stumbling blocks?

No, Thank you to my mentor that give me the right hints!

Report5 - 2012-06-29

What did I do this week?

Finish the RasterNumpy class, it's working now: 

>>> import pygrass
>>> elev = pygrass.RasterNumpy('elevation', 'PERMANENT')
>>> 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]
>>> # in this case RasterNumpy is an extention of the numpy class
>>> # therefore you may use all the fancy things of numpy.
>>> elev[:5, :3]
RasterNumpy([[ 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]], dtype=float32)
>>> el = elev < 144
>>> el[:5, :3]
RasterNumpy([[1, 1, 1],
       [1, 1, 1],
       [1, 1, 1],
       [0, 0, 0],
       [0, 0, 0]], dtype=int32)
>>> el.name == None
True
>>> # give a name to the new map
>>> el.name = 'new'
>>> el.close()
>>> el.remove()

Add the RasterNumpy class to the benchmarks. Start to study the C API to implement the categories.

What will I be working on next week?

Go on to develop categories support for the raster classes,

Did I meet with any stumbling blocks?

No

References

  1. ↑ http://numpy.scipy.org/
Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=GRASS_GSoC_2012_High_level_map_interaction&oldid=16013"