MapFish protocol modifications for vector simplification

protocol.py

#
# Copyright (C) 2009  Camptocamp
#
# This file is part of MapFish Server
#
# MapFish Server is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# MapFish Server is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with MapFish Server.  If not, see <http://www.gnu.org/licenses/>.
#

import logging
log = logging.getLogger(__name__)

from pylons.controllers.util import abort

from shapely.geometry import asShape
from shapely.geometry.point import Point
from shapely.geometry.polygon import Polygon

from sqlalchemy.sql import asc, desc, and_

from geoalchemy import WKBSpatialElement
from geoalchemy.functions import functions

from geojson import Feature, FeatureCollection, loads, GeoJSON

from mapfish.sqlalchemygeom import within_distance

import math

# Start GeoJSON Tiles
TOLERANCES = {
	17: 0.00001,
	16: 0.00002,
	15: 0.00004,
	14: 0.00008,
	13: 0.00016,
	12: 0.00032,
	11: 0.00064,
	10: 0.00128,
	9: 0.00256,
	8: 0.00512,
	7: 0.01024,
	6: 0.02048,
	5: 0.04096,
	4: 0.08182,
	3: 0.16384,
	2: 0.32768,
	1: 0.65536,
	0: 1.31072
}

def zoom2tolerance(zoom):
    if zoom < 0:
        return 1.31072
    if zoom > 17:
        return 0.00001
    return TOLERANCES[zoom]

def TileBounds(tx, ty, zoom):
	# Returns bounds of the given tile in EPSG:4326 coordinates
	
	ty = ((1 << zoom) - ty - 1)
	
	# Pixels to Meters
	# Converts pixel coordinates in given zoom level of pyramid to EPSG:900913
	originShift = 20037508.342789244
	# resolution = (meters/pixel) for given zoom level (measured at Equator)
	resolution = 156543.03392804062 / (2**zoom)
	minx_meters = ( tx * 256 ) * resolution - originShift
	miny_meters = ( ty * 256 ) * resolution - originShift
	maxx_meters = ( ( tx + 1 ) * 256 ) * resolution - originShift
	maxy_meters = ( ( ty + 1 ) * 256 ) * resolution - originShift

	# Meters to Lat/Lon
	minx_ll = (minx_meters / originShift) * 180.0
	miny_ll = (miny_meters / originShift) * 180.0
	maxx_ll = (maxx_meters / originShift) * 180.0
	maxy_ll = (maxy_meters / originShift) * 180.0
	
	miny_ll = 180 / math.pi * (2 * math.atan( math.exp( miny_ll * math.pi / 180.0)) - math.pi / 2.0)
	maxy_ll = 180 / math.pi * (2 * math.atan( math.exp( maxy_ll * math.pi / 180.0)) - math.pi / 2.0)
	
	
	return Polygon(((minx_ll, miny_ll), (minx_ll, maxy_ll), (maxx_ll, maxy_ll), (maxx_ll, miny_ll),(minx_ll, miny_ll)))
# End GeoJSON Tiles

def create_geom_filter(request, mapped_class, **kwargs):
    """Create MapFish geometry filter based on the request params. Either
    a box or within or geometry filter, depending on the request params.
    Additional named arguments are passed to the spatial filter."""

    tolerance = 0
    if 'tolerance' in request.params:
        tolerance = float(request.params['tolerance'])

    epsg = None
    if 'epsg' in request.params:
        epsg = int(request.params['epsg'])

    box = None
    if 'bbox' in request.params:
        box = request.params['bbox']

    geometry = None

    if box is not None:
        box = map(float, box.split(','))
        geometry = Polygon(((box[0], box[1]), (box[0], box[3]),
                            (box[2], box[3]), (box[2], box[1]),
                            (box[0], box[1])))
    elif 'lon' and 'lat' in request.params:
        geometry = Point(float(request.params['lon']),
                         float(request.params['lat']))
    elif 'geometry' in request.params:
        factory = lambda ob: GeoJSON.to_instance(ob)
        geometry = loads(request.params['geometry'], object_hook=factory)
        geometry = asShape(geometry)

    if geometry is None:
        return None

    geom_column = mapped_class.geometry_column()

    epsg = geom_column.type.srid if epsg is None else epsg
    if epsg != geom_column.type.srid:
        geom_column = functions.transform(geom_column, epsg)

    wkb_geometry = WKBSpatialElement(buffer(geometry.wkb), epsg)

    if 'additional_params' in kwargs:
        return within_distance(geom_column, wkb_geometry, tolerance,
                               kwargs['additional_params'])
    else:
        return within_distance(geom_column, wkb_geometry, tolerance)
 
def create_attr_filter(request, mapped_class):
    """Create an ``and_`` SQLAlchemy filter (a ClauseList object) based
    on the request params (``queryable``, ``eq``, ``ne``, ...)."""

    mapping = {
        'eq'   : '__eq__',
        'ne'   : '__ne__',
        'lt'   : '__lt__',
        'lte'  : '__le__',
        'gt'   : '__gt__',
        'gte'  : '__ge__',
        'like' : 'like',
        'ilike': 'ilike'
    }

    filters = []
    if 'queryable' in request.params:
        queryable = request.params['queryable'].split(',')
        for k in request.params:
            if len(request.params[k]) <= 0 or '__' not in k:
                continue

            col, op = k.split("__")
            if col not in queryable or op not in mapping.keys():
                continue

            column = getattr(mapped_class, col)
            f = getattr(column, mapping[op])(request.params[k])
            filters.append(f)

    return and_(*filters) if len(filters) > 0 else None

def create_default_filter(request, mapped_class, **kwargs):
    """ Create MapFish default filter based on the request params. Additional
    named arguments are passed to the spatial filter."""

    attr_filter = create_attr_filter(request, mapped_class)
    geom_filter = create_geom_filter(request, mapped_class, **kwargs)

    if geom_filter is None and attr_filter is None:
        return None

    return and_(geom_filter, attr_filter)

def asbool(val):
    # Convert the passed value to a boolean.
    if isinstance(val, str) or isinstance(val, unicode):
        low = val.lower()
        return low != 'false' and low != '0'
    else:
        return bool(val)

class Protocol(object):
    """ Protocol class.

      Session
          the SQLAlchemy session.

      mapped_class
          the class mapped to a database table in the ORM.

      readonly
          ``True`` if this protocol is read-only, ``False`` otherwise. If
          ``True``, the methods ``create()``, ``update()`` and  ``delete()``
          will set 403 as the response status and return right away.

      \**kwargs
          before_create
            a callback function called before a feature is inserted
            in the database table, the function receives the request,
            the feature read from the GeoJSON document sent in the
            request, and the database object to be updated. The
            latter is None if this is is an actual insertion.

          before_update
            a callback function called before a feature is updated
            in the database table, the function receives the request,
            the feature read from the GeoJSON document sent in the
            request, and the database object to be updated.

          before_delete
            a callback function called before a feature is deleted
            in the database table, the function receives the request
            and the database object about to be deleted.
    """

    def __init__(self, Session, mapped_class, readonly=False, **kwargs):

        self.Session = Session
        self.mapped_class = mapped_class
        self.readonly = readonly
        self.before_create = None
        if kwargs.has_key('before_create'):
            self.before_create = kwargs['before_create']
        self.before_update = None
        if kwargs.has_key('before_update'):
            self.before_update = kwargs['before_update']
        self.before_delete = None
        if kwargs.has_key('before_delete'):
            self.before_delete = kwargs['before_delete']

    def _filter_attrs(self, feature, request):
        """ Remove some attributes from the feature and set the geometry to None
            in the feature based ``attrs`` and the ``no_geom`` parameters. """
        if 'attrs' in request.params:
            attrs = request.params['attrs'].split(',')
            props = feature.properties
            new_props = {}
            for name in attrs:
                if name in props:
                    new_props[name] = props[name]
            feature.properties = new_props

        if asbool(request.params.get('no_geom', False)):
            feature.geometry=None
        return feature

    def _get_order_by(self, request):
        """ Return an SA order_by """
        column_name = None
        if 'sort' in request.params:
            column_name = request.params['sort']
        elif 'order_by' in request.params:
            column_name = request.params['order_by']

        if column_name and column_name in self.mapped_class.__table__.c:
            column = self.mapped_class.__table__.c[column_name]
            if 'dir' in request.params and request.params['dir'].upper() == 'DESC':
                return desc(column)
            else:
                return asc(column)
        else:
            return None

    def _query(self, request, filter=None, execute=True):
        """ Build a query based on the filter and the request params,
            and send the query to the database. """

        limit = None
        offset = None

        if 'maxfeatures' in request.params:
            limit = int(request.params['maxfeatures'])
        if 'limit' in request.params:
            limit = int(request.params['limit'])
        if 'offset' in request.params:
            offset = int(request.params['offset'])

        if filter is None:
            # create MapFish default filter
            filter = create_default_filter(request, self.mapped_class)

        query = self.Session.query(self.mapped_class).filter(filter)

        order_by = self._get_order_by(request)
        if order_by is not None:
            query = query.order_by(order_by)

        query = query.limit(limit).offset(offset)

        if execute:
            return query.all()
        else:
            return query

    def count(self, request, filter=None):
        """ Return the number of records matching the given filter. """
        if filter is None:
            filter = create_default_filter(request, self.mapped_class)
        return str(self.Session.query(self.mapped_class).filter(filter).count())

    def geojson(self, request, response, x=None, y=None, z=None):
        """ Return features in X/Y/Z tile bounds. """
        if not x and y and z:
            abort(404)
            
        x = int(x)
        y = int(y)
        z = int(z)    
        
        geom_column = self.mapped_class.geometry_column()
        #epsg = geom_column.type.srid
        geometry = TileBounds(x, y, z)
        wkb_geometry = WKBSpatialElement(buffer(geometry.wkb), 4326)

        geom_filter = within_distance(geom_column, wkb_geometry, 0)

        attr_filter = create_attr_filter(request, self.mapped_class)
        filter = and_(attr_filter, geom_filter)
        
        response.headers["Access-Control-Allow-Origin"] = "*"
        
        tolerance = zoom2tolerance(z)
        
        return self.read(request, filter, None, tolerance)

    def read(self, request, filter=None, id=None, tolerance=None):
        """ Build a query based on the filter or the idenfier, send the query
        to the database, and return a Feature or a FeatureCollection. """
        ret = None
        if id is not None:
            o = self.Session.query(self.mapped_class).get(id)
            if o is None:
                abort(404)
            ret = self._filter_attrs(o.toFeature(), request)
        else:
            objs = self._query(request, filter)
            if tolerance is None:
                ret = FeatureCollection(
                    [self._filter_attrs(o.toFeature(), request) \
                        for o in objs if o.geometry is not None])
            else:
                ret = FeatureCollection(
                    [self._filter_attrs(o.toFeature(tolerance=tolerance), request) \
                        for o in objs if o.geometry is not None])
        return ret

    def create(self, request, response, execute=True):
        """ Read the GeoJSON feature collection from the request body and
            create new objects in the database. """
        if self.readonly:
            abort(403)
        content = request.environ['wsgi.input'].read(int(request.environ['CONTENT_LENGTH']))
        factory = lambda ob: GeoJSON.to_instance(ob)
        collection = loads(content, object_hook=factory)
        if not isinstance(collection, FeatureCollection):
            abort(400)
        objects = []
        for feature in collection.features:
            create = False
            obj = None
            if isinstance(feature.id, int):
                obj = self.Session.query(self.mapped_class).get(feature.id)
            if self.before_create is not None:
                self.before_create(request, feature, obj)
            if obj is None:
                obj = self.mapped_class()
                create = True
            self.__copy_attributes(feature, obj)
            if create:
                self.Session.add(obj)
            objects.append(obj)
        if execute:
            self.Session.commit()
        response.status = 201
        if len(objects) > 0:
            features = [o.toFeature() for o in objects \
                            if o.geometry is not None]
            return FeatureCollection(features)
        return

    def update(self, request, response, id):
        """ Read the GeoJSON feature from the request body and update the
        corresponding object in the database. """
        if self.readonly:
            abort(403)
        obj = self.Session.query(self.mapped_class).get(id)
        if obj is None:
            abort(404)
        content = request.environ['wsgi.input'].read(int(request.environ['CONTENT_LENGTH']))
        factory = lambda ob: GeoJSON.to_instance(ob)
        feature = loads(content, object_hook=factory)
        if not isinstance(feature, Feature):
            abort(400)
        if self.before_update is not None:
            self.before_update(request, feature, obj)
        self.__copy_attributes(feature, obj)
        self.Session.commit()
        response.status = 201
        return obj.toFeature()

    def delete(self, request, response, id):
        """ Remove the targetted feature from the database """
        if self.readonly:
            abort(403)
        obj = self.Session.query(self.mapped_class).get(id)
        if obj is None:
            abort(404)
        if self.before_delete is not None:
            self.before_delete(request, obj)
        self.Session.delete(obj)
        self.Session.commit()
        response.status = 204
        return

    def __copy_attributes(self, json_feature, obj):
        """Updates the passed-in object with the values
        from the GeoJSON feature."""
        # create a Shapely geometry from GeoJSON and persist the geometry using WKB
        shape = asShape(json_feature.geometry)
        srid = self.mapped_class.geometry_column().type.srid
        obj.geometry = WKBSpatialElement(buffer(shape.wkb), srid=srid)
        # also store the Shapely geometry so that we can use it to return the geometry as GeoJSON
        obj.geometry.shape = shape
        for key in json_feature.properties:
            obj[key] = json_feature.properties[key]