diff --git a/conda/meta.yaml b/conda/meta.yaml
index fc4a40d..2d94fbe 100644
--- a/conda/meta.yaml
+++ b/conda/meta.yaml
@@ -1,4 +1,4 @@
-{% set version = "2.0.0" %}
+{% set version = "2.1.0" %}
 
 package:
   name: alinea.astk
diff --git a/src/alinea/astk/scene/__init__.py b/src/alinea/astk/scene/__init__.py
new file mode 100644
index 0000000..7ea52e6
--- /dev/null
+++ b/src/alinea/astk/scene/__init__.py
@@ -0,0 +1 @@
+"""A collection of simple scene objects that can be maped to/from PlanGL scenes"""
\ No newline at end of file
diff --git a/src/alinea/astk/scene/color_map.py b/src/alinea/astk/scene/color_map.py
new file mode 100644
index 0000000..6a01a5a
--- /dev/null
+++ b/src/alinea/astk/scene/color_map.py
@@ -0,0 +1,63 @@
+class ColorMap(object):
+    """A RGB color map, between 2 colors defined in HSV code
+    """
+
+    def __init__(self, minval=0., maxval=1.):
+        self.minval = float(minval)
+        self.maxval = float(maxval)
+
+    def color(self, normedU):
+        """
+
+        :param normedU: todo
+
+        """
+        inter = 1 / 5.
+        winter = int(normedU / inter)
+        a = (normedU % inter) / inter
+        b = 1 - a
+
+        if winter < 0:
+            col = (self.coul2, self.coul2, self.coul1)
+        elif winter == 0:
+            col = (self.coul2, self.coul2 * b + self.coul1 * a, self.coul1)
+        elif winter == 1:
+            col = (self.coul2, self.coul1, self.coul1 * b + self.coul2 * a)
+        elif winter == 2:
+            col = (self.coul2 * b + self.coul1 * a, self.coul1, self.coul2)
+        elif winter == 3:
+            col = (self.coul1, self.coul1 * b + self.coul2 * a, self.coul2)
+        elif winter > 3:
+            col = (self.coul1, self.coul2, self.coul2)
+        return (int(col[0]), int(col[1]), int(col[2]))
+
+    def greycolor(self, normedU):
+        """
+
+        :param normedU: todo
+        :returns: todo
+        """
+        return (int(255 * normedU), int(255 * normedU), int(255 * normedU))
+
+    def grey(self, u):
+        """
+        :param u:
+        :returns: todo
+        """
+        return self.greycolor(self.normU(u))
+
+    def normU(self, u):
+        """
+        :param u:
+        :returns: todo
+        """
+        if self.minval == self.maxval:
+            return 0.5
+        return (u - self.minval) / (self.maxval - self.minval)
+
+    def __call__(self, u, minval=0, maxval=1, coul1=80, coul2=20):
+        self.coul1 = coul1
+        self.coul2 = coul2
+        self.minval = float(minval)
+        self.maxval = float(maxval)
+        return self.color(self.normU(u))
diff --git a/src/alinea/astk/scene/display.py b/src/alinea/astk/scene/display.py
new file mode 100644
index 0000000..2ae73cf
--- /dev/null
+++ b/src/alinea/astk/scene/display.py
@@ -0,0 +1,65 @@
+import warnings
+from itertools import chain
+from math import isnan
+
+import alinea.astk.scene.pgl_scene as pgls
+from alinea.astk.scene.color_map import ColorMap
+
+
+def jet_colors(x, minval=None, maxval=None):
+    """ return jet colors associated to a vector of values"""
+    if minval is None:
+        minval = min(x)
+    if maxval is None:
+        maxval = max(x)
+    cmap = ColorMap()
+
+    return map(lambda v: cmap(v, minval, maxval, 250., 20.), x)
+
+
+def nan_to_zero(values):
+    return [0 if isnan(x) else x for x in values]
+
+
+def property_as_colors(a_property, minval=None, maxval=None, gamma=None):
+    """ transform a scalar property in a color property of same kind
+
+        property is a {shape_id: value} or {shape_id: list of values} dict
+    """
+
+    values = a_property.values()
+    if isinstance(values[0], list):
+        values = list(chain.from_iterable(values))
+    values = nan_to_zero(values)
+    if minval is None:
+        minval = min(values)
+    if maxval is None:
+        maxval = max(values)
+    if gamma is None:
+        gamma = 1
+    norm = 0.5
+    if minval != maxval:
+        norm = maxval - minval
+    values = map(lambda x: ((x - minval) / float(norm)) ** gamma, values)
+    colors = jet_colors(values, 0, 1)
+    color_property = {}
+    for k, v in a_property.iteritems():
+        if isinstance(v, list):
+            color_property[k] = []
+            for i in range(len(v)):
+                color_property[k].append(colors.pop(0))
+        else:
+            color_property[k] = colors.pop(0)
+
+    return color_property
+
+
+def display_property(scene_mesh, property, minval=None, maxval=None):
+    colors = property_as_colors(property, minval=minval, maxval=maxval)
+    if pgls.pgl_imported:
+        scene = pgls.from_scene_mesh(scene_mesh, colors)
+        return pgls.display(scene)
+    else:
+        warnings.warn('PlanGL not found, no display available!!')
+        return scene_mesh
+
diff --git a/src/alinea/astk/scene/pgl_scene.py b/src/alinea/astk/scene/pgl_scene.py
new file mode 100644
index 0000000..928ab67
--- /dev/null
+++ b/src/alinea/astk/scene/pgl_scene.py
@@ -0,0 +1,104 @@
+""" Interfaces for PlantGL scene / PlantGL display"""
+
+import numpy
+
+pgl_imported = False
+try:
+    import openalea.plantgl.all as pgl
+    pgl_imported = True
+except ImportError:
+    pgl = None
+
+
+def bbox(pgl_scene, scene_unit='m'):
+    """ Bounding box of a pgl scene"""
+    tesselator = pgl.Tesselator()
+    bbc = pgl.BBoxComputer(tesselator)
+    bbc.process(pgl_scene)
+    box = bbc.result
+    xmin, ymin, zmin = box.getXMin(), box.getYMin(), box.getZMin()
+    xmax, ymax, zmax = box.getXMax(), box.getYMax(), box.getZMax()
+    if scene_unit != 'm':
+        units = {'mm': 0.001, 'cm': 0.01, 'dm': 0.1, 'm': 1, 'dam': 10,
+                 'hm': 100,
+                 'km': 1000}
+        convert = units.get(scene_unit, 1)
+        xmin, ymin, zmin = numpy.array((xmin, ymin, zmin)) * convert
+        xmax, ymax, zmax = numpy.array((xmax, ymax, zmax)) * convert
+
+    return (xmin, ymin, zmin), (xmax, ymax, zmax)
+
+
+def shape_mesh(pgl_shape, tesselator=None):
+    if tesselator is None:
+        tesselator = pgl.Tesselator()
+    tesselator.process(pgl_shape)
+    tset = tesselator.result
+    return numpy.array(tset.pointList), numpy.array(tset.indexList)
+
+
+def as_scene_mesh(pgl_scene):
+    """ Transform a PlantGL scene / PlantGL shape dict to a scene_mesh"""
+    tesselator = pgl.Tesselator()
+
+    if isinstance(pgl_scene, pgl.Scene):
+        sm = {}
+
+        def _concat_mesh(mesh1,mesh2):
+            v1, f1 = mesh1
+            v2, f2 = mesh2
+            v = numpy.array(v1.tolist() + v2.tolist())
+            offset = len(v1)
+            f = numpy.array(f1.tolist() + [[i + offset, j + offset, k + offset] for i, j, k
+                               in f2.tolist()])
+            return v, f
+
+        for pid, pgl_objects in pgl_scene.todict().iteritems():
+            sm[pid] = reduce(_concat_mesh, [shape_mesh(pgl_object, tesselator) for pgl_object in
+                           pgl_objects])
+        return sm
+    elif isinstance(pgl_scene, dict):
+        return {sh_id: shape_mesh(sh,tesselator) for sh_id, sh in
+                pgl_scene.iteritems()}
+    else:
+        return pgl_scene
+
+
+def from_scene_mesh(scene_mesh, colors=None):
+    plant_color = (0, 180, 0)
+
+    if colors is None:
+        colors = {k: plant_color for k in scene_mesh}
+
+    scene = pgl.Scene()
+    for sh_id, mesh in scene_mesh.iteritems():
+        vertices, faces = mesh
+        if isinstance(colors[sh_id], tuple):
+            r, g, b = colors[sh_id]
+            color_list = [pgl.Color4(r, g, b, 0)] * len(faces)
+        else:
+            color_list = [pgl.Color4(r, g, b, 0) for r, g, b in colors[sh_id]]
+        shape = pgl.TriangleSet(vertices, faces)
+        shape.colorList = color_list
+        shape.colorPerVertex = False
+        shape.id = sh_id
+        scene += shape
+
+    return scene
+
+
+def display(scene):
+    """ display a scene"""
+    pgl.Viewer.display(scene)
+    return scene
+
+
+def unit_sphere_scene():
+    return pgl.Scene([pgl.Sphere()])
+
+
+def is_pgl_scene(scene):
+    if not pgl_imported:
+        return False
+    else:
+        return isinstance(scene, pgl.Scene)
\ No newline at end of file
diff --git a/src/alinea/astk/scene/surfacic_point_cloud.py b/src/alinea/astk/scene/surfacic_point_cloud.py
new file mode 100644
index 0000000..d6f158a
--- /dev/null
+++ b/src/alinea/astk/scene/surfacic_point_cloud.py
@@ -0,0 +1,221 @@
+import numpy
+import pandas
+from alinea.astk.scene.tagged_mesh import spherical, surface, normal, \
+    centroid, random_normals, equilateral, move_points
+
+
+class SurfacicPointCloud(object):
+    """A data structure for handling discrete spatial distribution of surface"""
+
+    units = {'mm': 0.001, 'cm': 0.01, 'dm': 0.1, 'm': 1, 'dam': 10, 'hm': 100,
+             'km': 1000}
+
+    def __init__(self, x, y, z, area, normals=None, shape_id=None,
+                 properties=None, scene_unit='m'):
+        """Instantiate a SurfacicPointCloud canopy from a list of points
+
+        Args:
+            x: (array-like) x-coordinate of surfacic elements
+            y: (array-like) y-coordinate of surfacic elements
+            z: (array-like) z-coordinate of surfacic elements
+            area: (array-like) areas of surfacic elements
+            normals: (array of 3-tuples): coordinates of vector normal to
+             surfacic point. If None, normals are randomly sampled.
+            shape_id: (array-like) vector of identifiers that allow tagging
+             points as elements of the same scene object. If None (default),
+            shape_id is set to its index in input arrays.
+            properties: (name: {sh_id: value}) optional additional data that
+             allow associating scalar property to points via their shape_id.
+            scene_unit: (string) the length unit used for inputs.
+
+        Returns:
+            a surfacic point cloud instance with all data converted to meter
+            (m)
+        """
+
+        x, y, z, area = map(lambda val: numpy.array(val, ndmin=1),
+                            (x, y, z, area))
+
+        try:
+            self.convert = self.units[scene_unit]
+        except KeyError:
+            print 'Warning, unit', scene_unit, 'not found, meter assumed'
+            self.convert = 1
+
+        if normals is None:
+            normals = random_normals(len(x))
+
+        if shape_id is None:
+            shape_id = range(len(x))
+
+        if properties is None:
+            properties = {}
+
+        assert len(x) == len(y) == len(z) == len(area) == len(normals) == len(
+            shape_id)
+
+        self.size = len(x)
+        self.x = x * self.convert
+        self.y = y * self.convert
+        self.z = z * self.convert
+        self.area = area * self.convert**2
+        self.point_id = range(len(x))
+        self.shape_id = shape_id
+        self.normals = normals
+        self.properties = properties
+
+    @staticmethod
+    def from_scene_mesh(scene_mesh, properties=None, scene_unit='m'):
+        """ Instantiation from a scene mesh dict
+
+        Args:
+            scene_mesh: a {shape_id: (vertices, faces)} dict
+            properties: (name:{shape_id: entity}) optional additional
+                named data associated to surfacic points.
+            properties: (name: {sh_id: value}) optional additional data that
+             allow associating scalar property to points via their shape_id.
+            scene_unit: (string) the length unit used for inputs.
+        """
+
+        areas, x, y, z, normals, shape_id = [[] for _ in range(6)]
+        for sh_id, (vertices, faces) in scene_mesh.iteritems():
+            areas += [surface(f, vertices) for f in faces]
+            xx, yy, zz = zip(*[centroid(f, vertices) for f in faces])
+            x += list(xx)
+            y += list(yy)
+            z += list(zz)
+            normals += [normal(f, vertices) for f in faces]
+            shape_id.extend([sh_id] * len(faces))
+        return SurfacicPointCloud(x=x, y=y, z=z, area=areas, normals=normals,
+                                  shape_id=shape_id,
+                                  properties=properties, scene_unit=scene_unit)
+
+    def as_data_frame(self, add_properties=True):
+        nx, ny, nz = zip(*self.normals)
+        df = pandas.DataFrame({'point_id': self.point_id, 'x': self.x, 'y': self.y,
+                               'z': self.z, 'shape_id': self.shape_id,
+                               'area': self.area,
+                               'norm_x': nx, 'norm_y': ny, 'norm_z': nz})
+
+        if add_properties and len(self.properties) > 0:
+            d = {'shape_id': self.shape_id}
+            for k, v in self.properties.iteritems():
+                d.update({k: [v[sh] for sh in self.shape_id]})
+            df = df.merge(pandas.DataFrame(d))
+
+        return df
+
+    @staticmethod
+    def from_data_frame(df):
+        d = df.to_dict('list')
+        cols = (
+            'x', 'y', 'z', 'area', 'shape_id', 'point_id', 'norm_x', 'norm_y',
+            'norm_z')
+
+        property_cols = [col for col in d if col not in cols]
+        properties = None
+        if len(property_cols) > 0:
+            dfp = df.loc[:,['shape_id'] + property_cols].groupby('shape_id').agg(
+                lambda x: x.iloc[0]).reset_index()
+            dfpd = dfp.to_dict('list')
+            properties = {k: dict(zip(dfpd['shape_id'], dfpd[k])) for k in
+                          property_cols}
+
+        normals = zip(d['norm_x'], d['norm_y'], d['norm_z'])
+
+        return SurfacicPointCloud(x=d['x'], y=d['y'], z=d['z'], area=d['area'],
+                                  shape_id=d['shape_id'], normals=normals,
+                                  properties=properties)
+
+    def shape_map(self):
+        return pandas.DataFrame(
+            {'point_id': self.point_id, 'shape_id': self.shape_id
+             })
+
+    def area_map(self):
+        return pandas.DataFrame(
+            {'point_id': self.point_id, 'area': self.area
+             })
+
+    def xyz_map(self):
+        return pandas.DataFrame(
+            {'point_id': self.point_id, 'x': self.x, 'y': self.y, 'z': self.z
+             })
+
+    def as_scene_mesh(self):
+        """ A simple mesh representation of the point cloud"""
+
+        scene = {}
+        df = self.as_data_frame(add_properties=False)
+        for sh_id, g in df.groupby('shape_id'):
+            vertices = []
+            faces = []
+            nf = 0
+            for ind, row in g.iterrows():
+                tri = equilateral(row.area)
+                rot = numpy.random.rand() * numpy.pi / 3
+                norm = row.norm_x, row.norm_y, row.norm_z
+                pos = row.x, row.y, row.z
+                pts = move_points(tri, pos, norm, rot)
+                vertices.extend(pts)
+                faces.append((3 * nf, 3 * nf + 1, 3 * nf + 2))
+                nf += 1
+            scene[sh_id] = vertices, faces
+        return scene
+
+    def as_triangle_scene(self):
+        """ A {id: [triangles,..] representation of the point cloud"""
+        scene = {}
+        df = self.as_data_frame(add_properties=False)
+        for sh_id, g in df.groupby('shape_id'):
+            triangles=[]
+            for ind, row in g.iterrows():
+                tri = equilateral(row.area)
+                rot = numpy.random.rand() * numpy.pi / 3
+                norm = row.norm_x, row.norm_y, row.norm_z
+                pos = row.x, row.y, row.z
+                pts = move_points(tri, pos, norm, rot)
+                triangles.append(pts)
+            scene[sh_id] = triangles
+        return scene
+
+    def save(self, path='surfacic_point_cloud.csv'):
+        """ Save a csv representation of the object
+        """
+        df = self.as_data_frame()
+        df.to_csv(path, index=False)
+        return path
+
+    @staticmethod
+    def load(path='surfacic_point_cloud.csv'):
+        df = pandas.read_csv(path)
+        return SurfacicPointCloud.from_data_frame(df)
+
+    def bbox(self):
+        return (self.x.min(), self.y.min(), self.z.min()), (self.x.max(), \
+               self.y.max(), self.z.max())
+
+    def inclinations(self):
+        """ inclinations angles of normals (degrees, positive)"""
+        theta, phi = zip(*map(spherical, self.normals))
+        inclin = numpy.degrees(theta)
+        inclin = numpy.where(inclin == 90, 90, inclin % 90)
+        df = pandas.DataFrame(
+            {'point_id': self.point_id, 'shape_id': self.shape_id, 'inclination': inclin})
+        return df
+
+    def subset(self, point_id=None, shape_id=None):
+        """return a surfacic point cloud that is a subset of the calling
+        object"""
+
+        df = self.as_data_frame()
+        if point_id is not None:
+            if not hasattr(point_id, '__len__'):
+                point_id = [point_id]
+            df = df.loc[df.point_id.isin(point_id), :]
+        if shape_id is not None:
+            if not hasattr(shape_id, '__len__'):
+                shape_id = [shape_id]
+            df = df.loc[df.shape_id.isin(shape_id), :]
+
+        return self.from_data_frame(df)
diff --git a/src/alinea/astk/scene/tagged_mesh.py b/src/alinea/astk/scene/tagged_mesh.py
new file mode 100644
index 0000000..d8b11a8
--- /dev/null
+++ b/src/alinea/astk/scene/tagged_mesh.py
@@ -0,0 +1,114 @@
+"""A simple python tagged scene mesh of the form {tag:(vertices, faces),...}"""
+import numpy
+
+
+def norm(vector):
+    x, y, z = vector
+    return numpy.sqrt(x ** 2 + y ** 2 + z ** 2)
+
+
+def normed(vector):
+    """ normalised coordinates of (0,point) vector
+    """
+    return numpy.array(vector) / norm(vector)
+
+
+def spherical(vector):
+    """ inclination (theta) and azimuth (phi) spherical angles"""
+    x, y, z = normed(vector)
+    theta = numpy.arccos(z)
+    phi = numpy.arctan2(y, x)
+    return theta, phi
+
+
+def cartesian(theta, phi):
+    """ cartesian coordinates of a unit vector with inclination theta and
+    azimuth phi"""
+    x = numpy.sin(theta) * numpy.cos(phi)
+    y = numpy.sin(theta) * numpy.sin(phi)
+    z = numpy.cos(theta)
+    return x, y, z
+
+
+def surface(face, vertices):
+    a, b, c = [numpy.array(vertices[i]) for i in face]
+    return norm(numpy.cross(b - a, c - a)) / 2.0
+
+
+def normal(face, vertices):
+    a, b, c = [numpy.array(vertices[i]) for i in face]
+    return normed(numpy.cross(b - a, c - a))
+
+
+def centroid(face, vertices):
+    points = [numpy.array(vertices[i]) for i in face]
+    x, y, z = zip(*points)
+    return numpy.mean(x), numpy.mean(y), numpy.mean(z)
+
+
+def random_normals(size=1):
+    theta = numpy.pi / 2 * numpy.random.random_sample(size)
+    phi = 2 * numpy.pi * numpy.random.random_sample(size)
+    return zip(*cartesian(theta, phi))
+
+
+def rotation_matrix(axis, theta):
+    """
+    Return the rotation matrix associated with counterclockwise rotation arround
+    the given axis by theta radians.
+    """
+    axis = numpy.asarray(axis)
+    axis = axis / norm(axis)
+    a = numpy.cos(theta / 2.0)
+    b, c, d = -axis * numpy.sin(theta / 2.0)
+    aa, bb, cc, dd = a * a, b * b, c * c, d * d
+    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
+    return numpy.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
+                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
+                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
+
+
+def rotate(points, rotation_mat):
+    return [numpy.dot(rotation_mat, p) for p in points]
+
+
+def move_points(xy_points, position = (0, 0, 0), orientation=(0, 0, 1),
+                rotation=None):
+    """ Move 2d xy points arround origin at 3D position with orientation
+
+        If rotation is not None, points are rotated arround z before move
+    """
+
+    theta, phi = spherical(orientation)
+    roty = rotation_matrix((0, 1, 0), theta)
+    rotz = rotation_matrix((0, 0, 1), phi)
+    x, y = zip(*xy_points)
+    z = [0] * len(x)
+    pts = zip(x, y, z)
+    if rotation:
+        rot = rotation_matrix((0, 0, 1), rotation)
+        pts = rotate(pts, rot)
+    px, py, pz = position
+
+    newx, newy, newz = zip(*rotate(rotate(pts, roty), rotz))
+    newx, newy, newz = px + numpy.array(newx), py + numpy.array(newy),\
+                       pz + numpy.array(newz)
+
+    return zip(newx, newy, newz)
+
+
+def equilateral(area=1):
+    """2d coordinates of an equilateral triangle centered on origin"""
+
+    edge = 2 * numpy.sqrt(area / numpy.sqrt(3))
+    height = numpy.sqrt(area * numpy.sqrt(3))
+
+    return (
+    (-edge / 2., -height / 3.), (edge / 2., -height / 3.), (0, 2 * height / 3.))
+
+
+def unit_square_mesh():
+    """A 2 triangle mesh representing ground unit area"""
+    vertices = ((0, 0, 0), (1, 0, 0), (1, 1, 0), (0, 1, 0))
+    faces = ((0, 1, 3), (1, 2, 3))
+    return {0: (vertices, faces)}
\ No newline at end of file
diff --git a/src/alinea/astk/version.py b/src/alinea/astk/version.py
index 837edc6..4cb5a34 100644
--- a/src/alinea/astk/version.py
+++ b/src/alinea/astk/version.py
@@ -1,8 +1,8 @@
 # {# pkglts, version
 #  -*- coding: utf-8 -*-
 
-major = 0
-minor = 3
+major = 2
+minor = 1
 post = 0
 
 __version__ = ".".join([str(s) for s in (major, minor, post)])
diff --git a/test/test_scene/__init__.py b/test/test_scene/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/test/test_scene/test_display.py b/test/test_scene/test_display.py
new file mode 100644
index 0000000..43e6593
--- /dev/null
+++ b/test/test_scene/test_display.py
@@ -0,0 +1,56 @@
+import numpy
+from alinea.astk.scene.display import jet_colors, property_as_colors, \
+    display_property
+
+
+def test_jet_colors():
+    colors = jet_colors(range(10))
+    expected = [(20, 20, 250),
+                (20, 147, 250),
+                (20, 250, 224),
+                (20, 250, 96),
+                (71, 250, 20),
+                (198, 250, 20),
+                (250, 173, 20),
+                (250, 45, 20),
+                (250, 20, 20),
+                (250, 20, 20)]
+    numpy.testing.assert_array_equal(colors, expected)
+
+    colors = jet_colors(range(10), minval=0, maxval=100)
+    expected = [(20, 20, 250),
+                (20, 31, 250),
+                (20, 43, 250),
+                (20, 54, 250),
+                (20, 66, 250),
+                (20, 77, 250),
+                (20, 89, 250),
+                (20, 100, 250),
+                (20, 111, 250),
+                (20, 123, 250)]
+
+    numpy.testing.assert_array_equal(colors, expected)
+
+
+def test_property_as_colors():
+    properties = {'a_property': {1: 9, 2: 0}, 'another': {1: [9], 2: [0, 5]}}
+    colors = property_as_colors(properties['a_property'])
+    assert colors[1] == (250, 20, 20)
+    assert colors[2] == (20, 20, 250)
+    colors = property_as_colors(properties['another'])
+    assert colors[2][0] == (20, 20, 250)
+    assert colors[2][1] == (198, 250, 20)
+
+    colors = property_as_colors(properties['a_property'], minval=0, maxval=100)
+    assert colors[1] == (20, 123, 250)
+    assert colors[2] == (20, 20, 250)
+
+
+def test_display_property():
+    faces = ((0, 1, 2), (0, 2, 3), (0, 1, 3))
+    vertices = ((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1))
+    s_mesh = {1: (vertices, (faces[1],)),
+              2: (vertices, [faces[i] for i in (0, 2)])}
+    a_property = {1: 0, 2: 9}
+    sc = display_property(s_mesh, a_property)
+
diff --git a/test/test_scene/test_geometry.py b/test/test_scene/test_geometry.py
new file mode 100644
index 0000000..54f469b
--- /dev/null
+++ b/test/test_scene/test_geometry.py
@@ -0,0 +1,36 @@
+import numpy
+from alinea.astk.scene.tagged_mesh import normed, spherical, cartesian, \
+    surface, normal, centroid, random_normals, move_points, unit_square_mesh
+
+
+def test_spherical():
+    v = normed((1, 0, 1))
+    theta, phi = spherical(v)
+    assert phi == 0
+    numpy.testing.assert_almost_equal(theta, numpy.pi / 4)
+    cart = cartesian(theta, phi)
+    numpy.testing.assert_almost_equal(cart, v)
+
+
+def test_triangle_math():
+    face = range(3)
+    vertices = ((0, 0, 0), (1, 0, 0), (0, 1, 0))
+
+    numpy.testing.assert_almost_equal(surface(face, vertices), 0.5)
+    numpy.testing.assert_almost_equal(normal(face, vertices), (0, 0, 1))
+    numpy.testing.assert_almost_equal(centroid(face, vertices), (1./3, 1./3, 0))
+
+
+def test_random_normals():
+    n = random_normals(1)
+    assert len(n) == 1
+    assert len(n[0]) == 3
+    n = random_normals(10)
+    assert len(n) == 10
+    assert len(n[0]) == 3
+
+
+def test_move_points():
+    vertices, faces = unit_square_mesh().values()[0]
+    x, y, z = zip(*vertices)
+    newpoints = move_points(zip(x, y))
\ No newline at end of file
diff --git a/test/test_scene/test_pgl_scene.py b/test/test_scene/test_pgl_scene.py
new file mode 100644
index 0000000..7bd833b
--- /dev/null
+++ b/test/test_scene/test_pgl_scene.py
@@ -0,0 +1,56 @@
+import alinea.astk.scene.pgl_scene as pgls
+import numpy
+
+if pgls.pgl_imported:
+
+    def test_bbox():
+        scene = pgls.unit_sphere_scene()
+        box = pgls.bbox(scene)
+        numpy.testing.assert_array_equal(box,
+                                         ((-0.5, -0.5, -0.5), (0.5, 0.5, 0.5)))
+        box = pgls.bbox(scene, scene_unit='dam')
+        numpy.testing.assert_array_equal(box, ((-5, -5, -5), (5, 5, 5)))
+
+
+    def test_shape_mesh():
+        shape = pgls.pgl.Sphere()
+        vertices, faces = pgls.shape_mesh(shape)
+        assert len(vertices) == 58
+        assert len(vertices[0]) == 3
+        assert len(faces) == 112
+        assert len(faces[0]) == 3
+
+    def test_as_scene_mesh():
+        sphere = pgls.pgl.Sphere()
+        scene = pgls.pgl.Scene([sphere, sphere])
+        sh_id = [sh.id for sh in scene]
+        s_mesh = pgls.as_scene_mesh(scene)
+        assert len(s_mesh) == 2
+        assert len(s_mesh[sh_id[1]]) == 2
+
+        geom = {'bud': sphere}
+        s_mesh = pgls.as_scene_mesh(geom)
+        assert 'bud' in s_mesh
+
+
+    def test_from_scene_mesh():
+        faces = ((0, 1, 2), (0, 2, 3), (0, 1, 3))
+        vertices = ((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1))
+        s_mesh = {1: (vertices, (faces[1],)),
+                  2: (vertices, [faces[i] for i in (0, 2)])}
+        # default color
+        scene = pgls.from_scene_mesh(s_mesh)
+        assert len(scene) == 2
+        assert scene[0].geometry.colorList[0].green == 180
+        # one color per shape
+        scene = pgls.from_scene_mesh(s_mesh,
+                                     colors={1: (0, 0, 0), 2: (1, 1, 1)})
+        assert scene[0].geometry.colorList[0].green == 0
+        assert scene[1].geometry.colorList[0].green == 1
+        # one color per triangle
+        scene = pgls.from_scene_mesh(s_mesh,
+                                     colors={1: (0, 0, 0),
+                                             2: [(1, 1, 1), (2, 2, 2)]})
+        assert scene[0].geometry.colorList[0].green == 0
+        assert scene[1].geometry.colorList[0].green == 1
+        assert scene[1].geometry.colorList[1].green == 2
diff --git a/test/test_scene/test_surfacic_point_cloud.py b/test/test_scene/test_surfacic_point_cloud.py
new file mode 100644
index 0000000..3d0290c
--- /dev/null
+++ b/test/test_scene/test_surfacic_point_cloud.py
@@ -0,0 +1,107 @@
+import os
+import tempfile
+
+import numpy
+from alinea.astk.scene.surfacic_point_cloud import SurfacicPointCloud
+
+
+def test_spc_instantiation():
+    spc = SurfacicPointCloud(0, 0, 0, 1)
+    for w in (
+    'x', 'y', 'z', 'area', 'shape_id', 'normals', 'properties', 'size'):
+        assert hasattr(spc, w)
+    spc = SurfacicPointCloud(100, 100, 100, 10000, scene_unit='cm')
+    assert spc.x == spc.y == spc.z == spc.area == 1  # m2
+
+    faces = (range(3),)
+    vertices = ((0, 0, 0), (1, 0, 0), (0, 1, 0))
+    sc = {'0': (vertices, faces)}
+    spc = SurfacicPointCloud.from_scene_mesh(sc)
+    assert spc.area == 0.5
+    numpy.testing.assert_array_equal(spc.normals, ((0, 0, 1),))
+
+
+def test_data_frame():
+    spc = SurfacicPointCloud(0, 0, 0, 1)
+    df = spc.as_data_frame()
+    assert df.shape == (1, 9)
+    spc.properties.update({'a_property': {0: 3}})
+    df = spc.as_data_frame()
+    assert df.shape == (1, 10)
+
+
+def test_as_scene_mesh():
+    spc = SurfacicPointCloud(0, 0, 0, 1)
+    sc = spc.as_scene_mesh()
+    assert 0 in sc
+    v, f, = sc[0]
+    assert len(f) == 1
+    assert len(v) == 3
+
+
+def test_as_triangle_scene():
+    spc = SurfacicPointCloud(0, 0, 0, 1)
+    sc = spc.as_triangle_scene()
+    assert 0 in sc
+    triangles = sc[0]
+    assert len(triangles) == 1
+    pts = triangles[0]
+    assert len(pts) == 3
+    assert len(pts[0]) == 3
+
+
+def test_serialisation():
+    spc = SurfacicPointCloud(0, 0, 0, 1)
+    try:
+        tmpdir = tempfile.mkdtemp()
+        path = os.path.join(tmpdir, 'test.csv')
+        spc.save(path)
+        spc2 = SurfacicPointCloud.load(path)
+        assert spc.x == spc2.x
+        numpy.testing.assert_almost_equal(spc.normals, spc2.normals)
+
+        spc.properties.update({'a_property': {0: 3}})
+        spc.save(path)
+        spc2 = SurfacicPointCloud.load(path)
+        assert 'a_property' in spc2.properties
+        for k in spc.properties['a_property']:
+            assert k in spc2.properties['a_property']
+    except Exception as e:
+        raise e
+    finally:
+        os.remove(path)
+        os.rmdir(tmpdir)
+
+
+def test_bbox():
+    faces = ((0, 1, 2), (0, 2, 3), (0, 1, 3))
+    vertices = ((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1))
+    sc = {'0': (vertices, faces)}
+    spc = SurfacicPointCloud.from_scene_mesh(sc)
+    expected = ((0.0, 0.0, 0.0), [1. / 3] * 3)
+    numpy.testing.assert_almost_equal(spc.bbox(), expected)
+
+
+def test_inclinations():
+    faces = ((0, 1, 2), (0, 2, 3), (0, 1, 3))
+    vertices = ((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1))
+    sc = {1: (vertices, (faces[1],)), 2: (vertices, [faces[i] for i in (0, 2)])}
+    spc = SurfacicPointCloud.from_scene_mesh(sc)
+    df = spc.inclinations()
+    inc = {g: x.to_dict('list')['inclination'] for g, x in
+     df.groupby('shape_id')}
+    numpy.testing.assert_array_equal(inc[1], [90.0])
+    numpy.testing.assert_array_equal(inc[2], [0.0, 90.0])
+
+
+def test_subset():
+    faces = ((0, 1, 2), (0, 2, 3), (0, 1, 3))
+    vertices = ((0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1))
+    sc = {1: (vertices, (faces[1],)), 2: (vertices, [faces[i] for i in (0, 2)])}
+    spc = SurfacicPointCloud.from_scene_mesh(sc)
+    sub = spc.subset(point_id=1)
+    assert len(sub.as_data_frame()==1)
+    sub = spc.subset(shape_id=2)
+    assert len(sub.as_data_frame()==2)
+
+