def get_height(self, points, edge=True):
        """
        Given a set of points, computes the z value for the parametric
        equation of the Polygons in the Surface.
        
        This method computes recursively the ``Polygon.get_height()``
        method for all the Polygons in the Surface, obtaining the z 
        value for the points according to the local Polygon they belong.

        The points outside the object will have a NaN value in the
        z column. If the inputed points has a third column the z values
        outside the Surface's domain will remain unchanged, the rest
        will be replaced.

        :param points: list of coordinates of the points to calculate.
        :type points: ndarray (shape=(N, 2 or 3))
        :param edge: If True, consider the points in the Polygon's edge
            inside the Polygon.
        :type edge: bool
        :returns: (x, y, z) arrays
        :rtype: ndarray (shape=(N, 3))
        """
        for poly in self:
            points = poly.get_height(points, edge=edge, full=True)
        for hole in self.holes:
            pip = hole.contains(points, edge=True)
            points[pip, 2] = np.nan
        return points