thuban 1.2.2-12build3 / usr / lib / thuban / Thuban / Model / layer.py

# Copyright (c) 2001, 2002, 2003, 2004, 2005, 2007 by Intevation GmbH
# Authors:
# Bernhard Herzog <bh@intevation.de>
# Jonathan Coles <jonathan@intevation.de>
# Silke Reimer <silke@intevation.de>
#
# This program is free software under the GPL (>=v2)
# Read the file COPYING coming with Thuban for details.

__version__ = "$Revision: 2804 $"

import os
import warnings

from wxproj import shape_centroid

from Thuban import _

from messages import LAYER_PROJECTION_CHANGED, LAYER_VISIBILITY_CHANGED, \
     LAYER_CHANGED, LAYER_SHAPESTORE_REPLACED, CLASS_CHANGED

import classification

from color import Transparent, Black
from base import TitledObject, Modifiable
from data import SHAPETYPE_POLYGON, SHAPETYPE_ARC, SHAPETYPE_POINT
from label import ALIGN_CENTER, ALIGN_TOP, ALIGN_BOTTOM, \
     ALIGN_LEFT, ALIGN_RIGHT

import resource

from color import Color

shapetype_names = {SHAPETYPE_POINT: "Point",
                   SHAPETYPE_ARC: "Arc",
                   SHAPETYPE_POLYGON: "Polygon"}

class BaseLayer(TitledObject, Modifiable):

    """Base class for the layers."""

    def __init__(self, title, visible = True, projection = None):
        """Initialize the layer.

        title -- the title
        visible -- boolean. If true the layer is visible.
        """
        TitledObject.__init__(self, title)
        Modifiable.__init__(self)
        self.visible = visible
        self.projection = projection

    def Visible(self):
        """Return true if layer is visible"""
        return self.visible

    def SetVisible(self, visible):
        """Set the layer's visibility."""
        self.visible = visible
        self.issue(LAYER_VISIBILITY_CHANGED, self)

    def HasClassification(self):
        """Determine if this layer supports classifications."""
        return False

    def HasShapes(self):
        """Determine if this layer supports shapes."""
        return False

    def GetProjection(self):
        """Return the layer's projection."""
        return self.projection

    def SetProjection(self, projection):
        """Set the layer's projection."""
        self.projection = projection
        self.changed(LAYER_PROJECTION_CHANGED, self)

    def Type(self):
        return "Unknown"

class Layer(BaseLayer):

    """Represent the information of one geodata file (currently a shapefile)

    All children of the layer have the same type.

    A layer has fill and stroke colors. Colors should be instances of
    Color. They can also be Transparent, indicating no fill or no stroke.

    The layer objects send the following events, all of which have the
    layer object as parameter:

        TITLE_CHANGED -- The title has changed.
        LAYER_PROJECTION_CHANGED -- the projection has changed.
    """

    def __init__(self, title, data, projection = None,
                 fill = Transparent, 
                 stroke = Black, 
                 lineWidth = 1, 
                 visible = True):
        """Initialize the layer.

        title -- the title
        data -- datastore object for the shape data shown by the layer
        projection -- the projection object. Its Inverse method is
               assumed to map the layer's coordinates to lat/long
               coordinates
        fill -- the fill color or Transparent if the shapes are 
                not filled
        stroke -- the stroke color or Transparent if the shapes 
                are not stroked
        visible -- boolean. If true the layer is visible.

        colors are expected to be instances of Color class
        """
        BaseLayer.__init__(self, title, 
                                 visible = visible,
                                 projection = projection)

        self.__classification = None
        self.store = None

        self.SetShapeStore(data)

        self.classification_column = None
        self.SetClassificationColumn(None)
        self.SetClassification(None)

        self.__classification.SetDefaultLineColor(stroke)
        self.__classification.SetDefaultLineWidth(lineWidth)
        self.__classification.SetDefaultFill(fill)

        self.UnsetModified()

    def SetShapeStore(self, store):
        # Set the classification to None if there is a classification
        # and the new shapestore doesn't have a table with a suitable
        # column, i.e one with the same name and type as before
        # FIXME: Maybe we should keep it the same if the type is
        # compatible enough such as FIELDTYPE_DOUBLE and FIELDTYPE_INT
        if self.__classification is not None:
            columnname = self.classification_column
            columntype = self.GetFieldType(columnname)
            table = store.Table()
            if (columnname is not None
                and (not table.HasColumn(columnname)
                     or table.Column(columnname).type != columntype)):
                self.SetClassification(None)

        self.store = store

        self.changed(LAYER_SHAPESTORE_REPLACED, self)

    def ShapeStore(self):
        return self.store

    def Destroy(self):
        BaseLayer.Destroy(self)
        if self.__classification is not None:
            self.__classification.Unsubscribe(CLASS_CHANGED,
                                              self._classification_changed)

    def BoundingBox(self):
        """Return the layer's bounding box in the intrinsic coordinate system.

        If the layer has no shapes, return None.
        """
        return self.store.BoundingBox()

    def LatLongBoundingBox(self):
        """Return the layer's bounding box in lat/long coordinates.

        Return None, if the layer doesn't contain any shapes.
        """
        bbox = self.BoundingBox()
        if bbox is not None and self.projection is not None:
            bbox = self.projection.InverseBBox(bbox)
        return bbox

    def Type(self):
        return self.ShapeType();

    def ShapesBoundingBox(self, shapes):
        """Return a bounding box in lat/long coordinates for the given
        list of shape ids.

        If shapes is None or empty, return None.
        """

        if shapes is None or len(shapes) == 0: return None

        xs = []
        ys = []

        for id in shapes:
            bbox = self.Shape(id).compute_bbox()
            if self.projection is not None:
                bbox = self.projection.InverseBBox(bbox)
            left, bottom, right, top = bbox
            xs.append(left); xs.append(right)
            ys.append(bottom); ys.append(top)

        return (min(xs), min(ys), max(xs), max(ys))


    def GetFieldType(self, fieldName):
        if self.store:
            table = self.store.Table()
            if table.HasColumn(fieldName):
                return table.Column(fieldName).type
        return None

    def HasShapes(self):
        return True

    def NumShapes(self):
        """Return the number of shapes in the layer"""
        return self.store.NumShapes()

    def ShapeType(self):
        """Return the type of the shapes in the layer.

        The return value is one of the SHAPETYPE_* constants defined in
        Thuban.Model.data.
        """
        return self.store.ShapeType()

    def Shape(self, index):
        """Return the shape with index index"""
        return self.store.Shape(index)

    def ShapesInRegion(self, bbox):
        """Return an iterable over the shapes that overlap the bounding box.

        The bbox parameter should be the bounding box as a tuple in the
        form (minx, miny, maxx, maxy) in unprojected coordinates.
        """
        if self.projection is not None:
            # Ensure that region lies within the layer's bounding box
            # Otherwise projection of the region would lead to incorrect
            # values.
            clipbbox = self.__mangle_bounding_box(bbox)
            bbox = self.projection.ForwardBBox(clipbbox)
        return self.store.ShapesInRegion(bbox)

    def GetClassificationColumn(self):
        return self.classification_column

    def SetClassificationColumn(self, column):
        """Set the column to classifiy on, or None. If column is not None
        and the column does not exist in the table, raise a ValueError.
        """
        if column:
            columnType = self.GetFieldType(column)
            if columnType is None:
                raise ValueError()
        changed = self.classification_column != column
        self.classification_column = column
        if changed:
            self.changed(LAYER_CHANGED, self)

    def HasClassification(self):
        return True

    def GetClassification(self):
        return self.__classification

    def SetClassification(self, clazz):
        """Set the classification used by this layer to 'clazz'

        If 'clazz' is None a default classification is created.

        This issues a LAYER_CHANGED event.
        """

        if self.__classification is not None:
            self.__classification.Unsubscribe(CLASS_CHANGED,
                                              self._classification_changed)

        if clazz is None:
            clazz = classification.Classification()

        self.__classification = clazz
        self.__classification.Subscribe(CLASS_CHANGED,
                                        self._classification_changed)

        self._classification_changed()

    def _classification_changed(self):
        """Called from the classification object when it has changed."""
        self.changed(LAYER_CHANGED, self)
 
    def TreeInfo(self):
        items = []

        items.append(_("Filename: %s") % self.ShapeStore().FileName())

        if self.Visible():
            items.append(_("Shown"))
        else:
            items.append(_("Hidden"))
        items.append(_("Shapes: %d") % self.NumShapes())

        bbox = self.LatLongBoundingBox()
        if bbox is not None:
            items.append(_("Extent (lat-lon): (%g, %g, %g, %g)") % tuple(bbox))
        else:
            items.append(_("Extent (lat-lon):"))
        items.append(_("Shapetype: %s") % shapetype_names[self.ShapeType()])

        if self.projection and len(self.projection.params) > 0:
            items.append((_("Projection"),
                        [str(param) for param in self.projection.params]))

        items.append(self.__classification)

        return (_("Layer '%s'") % self.Title(), items)

    def __mangle_bounding_box(self, bbox):
        # FIXME: This method doesn't make much sense.
        # See RT #2845 which effectively says:
        #
        # If this method, which was originally called ClipBoundingBox,
        # is supposed to do clipping it shouldn't return the parameter
        # unchanged when it lies completely outside of the bounding box.
        # It would be better to return None and return an empty list in
        # ShapesInRegion (the only caller) in that case.
        #
        # This method was introduced to fix a bug that IIRC had
        # something todo with projections and bounding boxes containing
        # NaN or INF when the parameter to ShapesInRegion covered the
        # entire earth or something similarly large).
        bminx, bminy, bmaxx, bmaxy = bbox
        lminx, lminy, lmaxx, lmaxy = self.LatLongBoundingBox()
        if bminx > lmaxx or bmaxx < lminx:
            left, right = bminx, bmaxx
        else:
            left = max(lminx, bminx)
            right = min(lmaxx, bmaxx)
        if bminy > lmaxy or bmaxy < lminy:
            bottom, top = bminy, bmaxy
        else:
            bottom = max(lminy, bminy)
            top = min(lmaxy, bmaxy)

        return (left, bottom, right, top)

    def GetLabelPosFromShape(self, cmap, shape_index):
        '''
        Return the label position parameters (x, y, halign, valign) from the
        shape object
        '''
        proj = cmap.projection
        if proj is not None:
            map_proj = proj
        else:
            map_proj = None
        proj = self.projection
        if proj is not None:
            layer_proj = proj
        else:
            layer_proj = None

        shapetype = self.ShapeType()
        if shapetype == SHAPETYPE_POLYGON:
            shapefile = self.ShapeStore().Shapefile().cobject()
            x, y = shape_centroid(shapefile, shape_index,
                                  map_proj, layer_proj, 1, 1, 0, 0)
            if map_proj is not None:
                x, y = map_proj.Inverse(x, y)
        else:
            shape = self.Shape(shape_index)
            if shapetype == SHAPETYPE_POINT:
                x, y = shape.Points()[0][0]
            else:
                # assume SHAPETYPE_ARC
                points = shape.Points()[0]
                x, y = points[len(points) / 2]
            if layer_proj is not None:
                x, y = layer_proj.Inverse(x, y)
        if shapetype == SHAPETYPE_POINT:
            halign = ALIGN_LEFT
            valign = ALIGN_CENTER
        elif shapetype == SHAPETYPE_POLYGON:
            halign = ALIGN_CENTER
            valign = ALIGN_CENTER
        elif shapetype == SHAPETYPE_ARC:
            halign = ALIGN_LEFT
            valign = ALIGN_CENTER
        
        return (x, y, halign, valign)



if resource.has_gdal_support():
    from osgeo import gdal
    from osgeo.gdalconst import GA_ReadOnly

class RasterLayer(BaseLayer):

    MASK_NONE  = 0
    MASK_BIT   = 1
    MASK_ALPHA = 2

    def __init__(self, title, filename, projection = None, 
                 visible = True, opacity = 1, masktype = MASK_ALPHA):
        """Initialize the Raster Layer.

        title -- title for the layer.

        filename -- file name of the source image.

        projection -- Projection object describing the projection which
                      the source image is in.

        visible -- True is the layer should initially be visible.

        Throws IOError if the filename is invalid or points to a file that
        is not in a format GDAL can use.
        """

        BaseLayer.__init__(self, title, visible = visible)

        self.projection = projection
        self.filename = os.path.abspath(filename)

        self.bbox = -1

        self.mask_type = masktype
        self.opacity = opacity

        self.image_info = None

        if resource.has_gdal_support():
            #
            # temporarily open the file so that GDAL can test if it's valid.
            #
            dataset = gdal.Open(self.filename, GA_ReadOnly)

            if dataset is None:
                raise IOError()

            #
            # while we have the file, extract some basic information 
            # that we can display later
            #
            self.image_info = {}

            self.image_info["nBands"] = dataset.RasterCount
            self.image_info["Size"] = (dataset.RasterXSize, dataset.RasterYSize)
            self.image_info["Driver"] = dataset.GetDriver().ShortName

            # store some information about the individual bands
            # [min_value, max_value]
            a = self.image_info["BandData"] = []

            for i in range(1, dataset.RasterCount+1):
                band = dataset.GetRasterBand(i)
                a.append(band.ComputeRasterMinMax())

        self.UnsetModified()

    def BoundingBox(self):
        """Return the layer's bounding box in the intrinsic coordinate system.

        If the there is no support for images, or the file cannot
        be read, or there is no geographics information available, return None.
        """
        if not resource.has_gdal_support():
            return None

        if self.bbox == -1:
            dataset = gdal.Open(self.filename, GA_ReadOnly)
            if dataset is None:
                self.bbox = None
            else:
                geotransform = dataset.GetGeoTransform()
                if geotransform is None:
                    return None

                x = 0
                y = dataset.RasterYSize
                left = geotransform[0] +        \
                       geotransform[1] * x +    \
                       geotransform[2] * y

                bottom = geotransform[3] +      \
                         geotransform[4] * x +  \
                         geotransform[5] * y

                x = dataset.RasterXSize
                y = 0
                right = geotransform[0] +       \
                        geotransform[1] * x +   \
                        geotransform[2] * y

                top = geotransform[3] +         \
                      geotransform[4] * x +     \
                      geotransform[5] * y

                self.bbox = (left, bottom, right, top)

        return self.bbox

    def LatLongBoundingBox(self):
        bbox = self.BoundingBox()
        if bbox is None:
            return None

        if self.projection is not None:
            bbox = self.projection.InverseBBox(bbox)

        return bbox

    def Type(self):
        return "Image"

    def GetImageFilename(self):
        return self.filename

    def MaskType(self):
        """Return True if the mask should be used when rendering the layer."""
        return self.mask_type

    def SetMaskType(self, type):
        """Set the type of mask to use.

        type can be one of MASK_NONE, MASK_BIT, MASK_ALPHA

        If the state changes, a LAYER_CHANGED message is sent.
        """
        if type not in (self.MASK_NONE, self.MASK_BIT, self.MASK_ALPHA):
            raise ValueError("type is invalid")

        if type != self.mask_type: 
            self.mask_type = type
            self.changed(LAYER_CHANGED, self)

    def Opacity(self):
        """Return the level of opacity used in alpha blending.
        """
        return self.opacity

    def SetOpacity(self, op):
        """Set the level of alpha opacity. 

        0 <= op <= 1.

        The layer is fully opaque when op = 1.
        """
        if not (0 <= op <= 1):
            raise ValueError("op out of range")

        if op != self.opacity: 
            self.opacity = op
            self.changed(LAYER_CHANGED, self)

    def ImageInfo(self):
        return self.image_info

    def TreeInfo(self):
        items = []

        items.append(_("Filename: %s") % self.GetImageFilename())

        if self.Visible():
            items.append(_("Shown"))
        else:
            items.append(_("Hidden"))

        bbox = self.LatLongBoundingBox()
        if bbox is not None:
            items.append(_("Extent (lat-lon): (%g, %g, %g, %g)") % bbox)
        else:
            items.append(_("Extent (lat-lon):"))

        if self.projection and len(self.projection.params) > 0:
            items.append((_("Projection"),
                        [str(param) for param in self.projection.params]))

        return (_("Layer '%s'") % self.Title(), items)