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

See GRASS SoC Ideas 2012/High level map interaction/Vector for more info.

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	√
June 29	Add RasterNumpy class	√
July 6	Add categories to Rasters classes	√
July 13	Implement module wrapper class or function generator	√
July 13	Mid-term evaluations deadline	√
July 20	Continue the module wrapper part	√
July 27	Add basic vector features	√
August 3	Add complex vector features	√
August 10	Add support for vector attributes	√
August 13	Suggested 'pencils down' date. Take a week to scrub code, write tests, improve documentation, etc.	√
August 17	Add support to the Vector head structure	√
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

Report6 - 2012-07-06

What did I do this week?

Add categories to the Raster classes:

All the raster classes support raster categories and share commons methods to modify the raster category. It is possible to check if the map has or not the categories with the has_cats method.

>>> elev.has_cats()
False

Opening a map that has category, for example the “landcove_1m” raster map from the North Carolina mapset. The has_cats method return True.

>>> land = pygrass.RasterRow('landcover_1m')
>>> land.has_cats()
True

Get and set the categories title, with:

>>> land.cats_title
'Rural area: Landcover'
>>> land.cats_title = 'Rural area: Landcover2'
>>> land.cats_title
'Rural area: Landcover2'
>>> land.cats_title = 'Rural area: Landcover'

Get the number of categories of the map with:

>>> land.num_cats()
11

See all the categories with:

>>> land.cats
[('pond', 1, None),
 ('forest', 2, None),
 ('developed', 3, None),
 ('bare', 4, None),
 ('paved road', 5, None),
 ('dirt road', 6, None),
 ('vineyard', 7, None),
 ('agriculture', 8, None),
 ('wetland', 9, None),
 ('bare ground path', 10, None),
 ('grass', 11, None)]

Access to single category, using Rast_get_ith_cat(), with:

>>> land.cats[0]
('pond', 1, None)
>>> land.cats['pond']
('pond', 1, None)
>>> land.get_cat(0)
('pond', 1, None)
>>> land.get_cat('pond')
('pond', 1, None)

Add new or change existing categories:

>>> land.set_cat('label', 1)
>>> land.get_cat('label')
('label', 1, None)
>>> land.set_cat('pond', 1, 1)
Sort categories, with:

>>> land.sort_cats() Copy categories from another raster map with:

>>> land.copy_cats(elev)

Read and Write:

>>> land.read_cats()
>>> #land.write_cats()

Get a Category object or set from a Category object:

>>> cats = land.get_cats()
>>> land.set_cats(cats)

Export and import from a file:

>>> land.write_cats_rules('land_rules.csv', ';')
>>> land.read_cats_rules('land_rules.csv', ';')

Add a new History class:

Instantiate a new object, with:

>>> hist = pygrass.History()

read the history of a map:

>>> hist.read('aspect')

Then it is possible to read and change the history attributes:

>>> hist.creator
'helena'
>>> hist.creator = 'pietro'
>>> hist.src1
'raster elevation file elev_ned10m'
>>> hist.src2
''
>>> hist.keyword
'generated by r.slope.aspect'
>>> hist.date
datetime.datetime(2006, 11, 7, 1, 11, 23)
>>> hist.mapset
'PERMANENT'
>>> hist.maptype
'raster'
>>> hist.title
'asp_ned10m'

What will I be working on next week?

Next week I will work on a function factory to wrap all the Grass modules and see and interact with module as python functions.

I would like something like:

>>> import pygrass.modules.raster as r
>>> elev, slope, aspect, dxx = r.slope_aspect(elevation='elevation', 
...                                           slope='slope',
...                                           aspect='aspect', 
...                                           dxx='dxx', overwrite = True)

or like:

>>> import pygrass
>>> import pygrass.modules.raster as r
>>> elev = pygrass.RasterRow('elevation')
>>> slope = pygrass.RasterRowIO('slope')
>>> aspect = pygrass.RasterSegment('aspect')
>>> dxx = pygrass.RasterNumpy('dxx')
>>> r.slope_aspect(elevation=elev, slope=slope aspect=aspect, dxx=dxx, overwrite = True)

The behavior and the implementation need to be discuss with my mentor. This is just an idea.

Did I meet with any stumbling blocks?

No

Report7 - 2012-07-13

What did I do this week?

Start to work on the pygrass modules class.

What will I be working on next week?

Continue with the the pygrass modules class.

Did I meet with any stumbling blocks?

No, I wrote a first implementation of the pygrass modules, but after discussion with my mentor I have started to rewrite the class.

Report8 - 2012-07-20

What did I do this week?

Finish the pygrass modules class and start to study the GRASS vector API. Yesterday and today I had some troubles with my OS installation and hard disk.

The Module class read the xml module description and instantiate a Module class, that is a callable object, user can change parameters, and the object check if the parameter are or not valid.

>>> from pygrass.modules import Module
>>> slope_aspect = Module("r.slope.aspect")
>>> slope_aspect(elevation='elevation', slope='slp',  aspect='asp',
...              format='percent', overwrite=True)

or it is possible to run all in one line, similar to run_command, with:

>>> slope_aspect = Module("r.slope.aspect", elevation='elevation',
...                        slope='slp',  aspect='asp',
...                        format='percent', overwrite=True)

The Module class allow user to run the command later, with:

>>> slope_aspect = Module("r.slope.aspect", elevation='elevation',
...                        slope='slp',  aspect='asp',
...                        format='percent', overwrite=True, run_=False)
>>> slope_aspect()

The class allow to manage the process, and read stdout and stderr, with:

>>> slope_aspect = Module('r.slope.aspect')
>>> slope_aspect(elevation='elevation', slope='slp', aspect='asp', overwrite=True, finish_=False)
>>> slope_aspect.popen.wait() # *.kill(), *.terminate()
0
>>> out, err = slope_aspect.popen.communicate()

for more info on the Module class look at: http://pygrass.readthedocs.org/en/latest/modules.html

What will I be working on next week?

Work on the Vector interface.

Did I meet with any stumbling blocks?

Not really, I need to better understand of the GRASS vector API.

Report9 - 2012-07-27

What did I do this week?

I've worked on vector API, in particular on geometry feature.

Point

• class Point(x, y, z=None, is2D=True)

>>> pnt = Point(0, 0)
>>> pnt.x
0.0
>>> pnt.y
0.0
>>> pnt.z
>>> pnt.is2D
True
>>> pnt.coords()
(0.0, 0.0)
>>> pnt.z = 0
>>> pnt.is2D
False
>>> pnt
Point(0.000000, 0.000000, 0.000000)
>>> print pnt
POINT(0.000000, 0.000000, 0.000000)

• buffer(dist=None, dist_x=None, dist_y=None, angle=0, round_=True, tol=0.1)

Return an Area object using the Vect_point_buffer2 C function. Creates buffer around the point (px, py).

• coords()

Return a tuple with the point coordinates.

>>> pnt = Point(0, 0)
>>> pnt.coords()
(0.0, 0.0)

If the point is 2D return a x, y tuple.

• distance(pnt)

Calculate distance of 2 points, using the Vect_points_distance C function, If one of the point have z == None, return the 2D distance.

>>> pnt0 = Point(0, 0, 0)
>>> pnt1 = Point(1, 0)
>>> pnt0.distance(pnt1)
1.0
>>> pnt1.z = 1
>>> pnt1
Point(1.000000, 0.000000, 1.000000)
>>> pnt0.distance(pnt1)
1.4142135623730951

The distance method require a :class:Point or a tuple with the coordinates.

• get_wkb()

Return a “well know binary” (WKB) geometry buffer

• get_wkt()

Return a “well know text” (WKT) geometry string.

>>> pnt = Point(0, 0)
>>> pnt.get_wkt()
'POINT(0.000000, 0.000000)'

Only POINT (2/3D) are supported, POINTM and POINT with: XYZM are not supported yet.

Line

class Line(points=None, mapinfo=None, lineid=None, field=None, is2D=True) Instantiate a new Line with a list of tuple, or with a list of Point.

>>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
>>> line                               
Line([Point(0.000000, 0.000000),
      Point(1.000000, 1.000000),
      Point(2.000000, 0.000000),
      Point(1.000000, -1.000000)])

• append(pnt)

Appends one point to the end of a line, using the Vect_append_point C function.

>>> line = Line()
>>> line.append((10, 100))
>>> line
Line([Point(10.000000, 100.000000)])
>>> line.append((20, 200))
>>> line
Line([Point(10.000000, 100.000000), Point(20.000000, 200.000000)])

Like python list.

• bbox()

Return the bounding box of the line, using Vect_line_box C function.

>>> line = Line([(0, 0), (0, 1), (2, 1), (2, 0)])
>>> bbox = line.bbox()
>>> bbox.north
1.0
>>> bbox.south
0.0
>>> bbox.east
2.0
>>> bbox.west
0.0
>>> bbox.top
0.0
>>> bbox.bottom
0.0

It is possible to access to the C struct of the object with: bbox.c_bbox

• buffer(dist=None, dist_x=None, dist_y=None, angle=0, round_=True, tol=0.1)

Return the buffer area around the line, using the Vect_line_buffer2 C function. Not implemented yet.

• delete(indx)

Remove the point in the index position.

>>> line = Line([(0, 0), (1, 1), (2, 2)])
>>> line.delete(-1)
>>> line
Line([Point(0.000000, 0.000000), Point(1.000000, 1.000000)])

• distance(pnt)

Return a tuple with:

 - the closest point on the line,
 - the distance between these two points,
 - distance of point from segment beginning
 - distance of point from line

The distance is compute using the Vect_line_distance C function.

• extend(line, forward=True)

Appends points to the end of a line.

It is possible to extend a line, give a list of points, or directly with a line_pnts struct.

If forward is True the line is extend forward otherwise is extend backward. The method use the Vect_append_points C function.

>>> line = Line([(0, 0), (1, 1)])
>>> line.extend( Line([(2, 2), (3, 3)]) )
>>> line                           
Line([Point(0.000000, 0.000000),
      Point(1.000000, 1.000000),
      Point(2.000000, 2.000000),
      Point(3.000000, 3.000000)])

Like python list, it is possible to extend a line, with another line or with a list of points.

• from_wkt(wkt)

Read a WKT string.

>>> line = Line()
>>> line.from_wkt("LINESTRING(0 0,1 1,1 2)")
>>> line                           
Line([Point(0.000000, 0.000000),
      Point(1.000000, 1.000000),
      Point(1.000000, 2.000000)])

• get_first_cat()

Fetches FIRST category number for given vector line and field, using the Vect_get_line_cat C function.

• get_pnt(distance, angle=0, slope=0)

Return a Point object on line in the specified distance, using the Vect_point_on_line C function. Raise a ValueError If the distance exceed the Line length.

>>> line = Line([(0, 0), (1, 1)])
>>> line.get_pnt(5)                           
Traceback (most recent call last):
    ...
ValueError: The distance exceed the lenght of the line, that is: 1.414214
>>> line.get_pnt(1)
Point(0.707107, 0.707107)

• get_wkt()

Return a Well Known Text string of the line.

>>> line = Line([(0, 0), (1, 1), (1, 2)])
>>> line.get_wkt()
'LINESTRING(0.000000 0.000000, 1.000000 1.000000, 1.000000 2.000000)'


* insert(indx, pnt)
Insert new point at index position and move all old points at that position and above up, using Vect_line_insert_point C function.
<source lang="python">
>>> line = Line([(0, 0), (1, 1)])
>>> line.insert(0, Point(1.000000, -1.000000) )
>>> line                           
Line([Point(1.000000, -1.000000),
      Point(0.000000, 0.000000),
      Point(1.000000, 1.000000)])

• length()

Calculate line length, 3D-length in case of 3D vector line, using Vect_line_length C function.

>>> line = Line([(0, 0), (1, 1), (0, 1)])
>>> line.length()
2.414213562373095

• length_geodesic()

Calculate line length, usig Vect_line_geodesic_length C function.

>>> line = Line([(0, 0), (1, 1), (0, 1)])
>>> line.length_geodesic()
2.414213562373095

• pop(indx)

Return the point in the index position and remove from the Line.

>>> line = Line([(0, 0), (1, 1), (2, 2)])
>>> midle_pnt = line.pop(1)
>>> midle_pnt
Point(1.000000, 1.000000)
>>> line
Line([Point(0.000000, 0.000000), Point(2.000000, 2.000000)])

• prune()

Remove duplicate points, i.e. zero length segments, using Vect_line_prune C function.

>>> line = Line([(0, 0), (1, 1), (1, 1), (2, 2)])
>>> line.prune()
>>> line                           
Line([Point(0.000000, 0.000000),
      Point(1.000000, 1.000000),
      Point(2.000000, 2.000000)])

• prune_thresh(threshold)

Remove points in threshold, using the Vect_line_prune_thresh C funtion.

• remove(pnt)

Delete point at given index and move all points above down, using Vect_line_delete_point C function.

>>> line = Line([(0, 0), (1, 1), (2, 2)])
>>> line.remove((2, 2))
>>> line[-1]
Point(1.000000, 1.000000)


* reset()
Reset line, using Vect_reset_line C function.
<source lang="python">
>>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
>>> len(line)
4
>>> line.reset()
>>> len(line)
0
>>> line
Line([])

• reverse()

Reverse the order of vertices, using Vect_line_reverse C function.

>>> line = Line([(0, 0), (1, 1), (2, 2)])
>>> line.reverse()
>>> line                           
Line([Point(2.000000, 2.000000),
      Point(1.000000, 1.000000),
      Point(0.000000, 0.000000)])

• segment(start, end)

Create line segment. using the Vect_line_segment C function. </source>

• toarray()

Return an array of coordinates.

>>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
>>> line.toarray()                 
array([[ 0.,  0.],
       [ 1.,  1.],
       [ 2.,  0.],
       [ 1., -1.]])

• tolist()

Return a list of tuple.

>>> line = Line([(0, 0), (1, 1), (2, 0), (1, -1)])
>>> line.tolist()
[(0.0, 0.0), (1.0, 1.0), (2.0, 0.0), (1.0, -1.0)]

What will I be working on next week?

Implement the other vector features.

Did I meet with any stumbling blocks?

Not really, but I need more time to assimilate the ideas of the GRASS vector API.

Report10 - 2012-08-03

What did I do this week?

Go on with the vector API, add more geometry feature like border, isle, and area, clean the code and study vector documentation and functions. This week I did an university exam, the last one for this summer.

What will I be working on next week?

Start to test the new vector classes add test and documentation

Did I meet with any stumbling blocks?

As the last week I'm studying more and more the vector API of GRASS, I would like to have more time... :-)

Report11 - 2012-08-10

What did I do this week?

This week I did a lot of works:

• move the documentation from ReadTheDocs to: http://www.ing.unitn.it/~zambelli/projects/pygrass in this way I'm able to use the sphinx automodule to document the API.
• Add more examples and documentation on vector geometry features: http://www.ing.unitn.it/~zambelli/projects/pygrass/vector.html
• Add a basic support for the vector attributes: http://www.ing.unitn.it/~zambelli/projects/pygrass/attributes.html

To run the above examples you need to use the North Carolina mapset, copy and build the topology, and copy the table attributes from dbf to sqlite and/or postgresql.

What will I be working on next week?

I will add the write mode to the Vector objects.

Did I meet with any stumbling blocks?

None

Report12 - 2012-08-17

What did I do this week?

• Re-organize the code structure;
• Add more documentation for each class, and when possible add docstring http://www.ing.unitn.it/~zambelli/projects/pygrass/;
• Add README, AUTHOR, NEWS files to the project root;
• Add the setup.py file to the project;
• Split the Vector class into Vector (without topology) and in VectTopo (with the topology support);
• Add vector head information.

What will I be working on next week?

Going on cleaning and fixing the code.

Did I meet with any stumbling blocks?

None

Download

Source code:

• git clone https://code.google.com/p/pygrass/
• download it (will be integrated in GRASS 7) - Integrated as of Nov. 2012.