passerelle/passerelle/apps/opengis/models.py

337 lines
14 KiB
Python

# passerelle - uniform access to multiple data sources and services
# Copyright (C) 2017 Entr'ouvert
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program 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 Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
import math
import xml.etree.ElementTree as ET
from HTMLParser import HTMLParser
import six
import pyproj
from django.core.cache import cache
from django.db import models
from django.http import HttpResponse
from django.utils.text import slugify
from django.utils.translation import ugettext_lazy as _
from passerelle.base.models import BaseResource
from passerelle.utils.api import endpoint
from passerelle.utils.jsonresponse import APIError
def build_dict_from_xml(elem):
d = {}
for child in elem.find('.'):
if child.tag.startswith('{'):
continue
attribute_name = slugify(child.tag).replace('-', '_')
if child.text and child.text.strip():
d[attribute_name] = child.text.strip()
else:
d[attribute_name] = build_dict_from_xml(child)
return d
PROJECTIONS = (
('EPSG:2154', _('EPSG:2154 (Lambert-93)')),
('EPSG:3857', _('EPSG:3857 (WGS 84 / Pseudo-Mercator)')),
('EPSG:3945', _('EPSG:3945 (CC45)')),
('EPSG:4326', _('EPSG:4326 (WGS 84)')),
)
class OpenGIS(BaseResource):
category = _('Geographic information system')
wms_service_url = models.URLField(_('Web Map Service (WMS) URL'), max_length=256, blank=True)
wfs_service_url = models.URLField(_('Web Feature Service (WFS) URL'), max_length=256, blank=True)
query_layer = models.CharField(_('Query Layer'), max_length=256)
projection = models.CharField(_('GIS projection'), choices=PROJECTIONS,
default='EPSG:3857', max_length=16)
search_radius = models.IntegerField(_('Radius for point search'), default=5)
attributes_mapping = (
('road', ('road', 'road_name', 'street', 'street_name', 'voie', 'nom_voie', 'rue')),
('city', ('city', 'city_name', 'town', 'town_name', 'commune', 'nom_commune', 'ville', 'nom_ville')),
('house_number', ('house_number', 'number', 'numero', 'numero_voie', 'numero_rue')),
('postcode', ('postcode', 'postalCode', 'zipcode', 'codepostal', 'cp', 'code_postal', 'code_post')),
('country', ('country', 'country_name', 'pays', 'nom_pays'))
)
class Meta:
verbose_name = _('OpenGIS')
def get_service_version(self, service_type, service_url, renew=False):
if not service_url:
raise APIError('no %s URL declared' % service_type)
cache_key = 'opengis-%s-%s-version' % (service_type, self.id)
if not renew:
service_version = cache.get(cache_key)
if service_version:
return service_version
response = self.requests.get(service_url, params={'request': 'GetCapabilities'})
element = ET.fromstring(response.content)
service_version = element.attrib.get('version')
# cache version number for an hour
cache.set(cache_key, service_version, 3600)
return service_version
def get_wms_service_version(self, renew=False):
return self.get_service_version('wms', self.wms_service_url, renew=renew)
def get_wfs_service_version(self, renew=False):
return self.get_service_version('wfs', self.wfs_service_url, renew=renew)
def check_status(self):
if self.wms_service_url:
response = self.requests.get(
self.wms_service_url,
params={
'service': 'WMS',
'request': 'GetCapabilities'
})
response.raise_for_status()
if self.wfs_service_url:
response = self.requests.get(
self.wfs_service_url,
params={
'service': 'WFS',
'request': 'GetCapabilities'
})
response.raise_for_status()
@endpoint(perm='can_access', description='Get features',
parameters={
'type_names': {
'description': _('Type of feature to query'),
'example_value': 'feature'
},
'property_name': {
'description': _('Property to list'),
'example_value': 'nom_commune'
},
'cql_filter': {
'description': _('CQL filter applied to the query'),
'example_value': 'commune=\'Paris\''
},
'filter_property_name': {
'description': _('Property by which to filter'),
'example_value': 'voie'
},
'q': {
'description': _('Filter value'),
'example_value': 'rue du chateau'
},
'case-insensitive': {
'description': _('Enables case-insensitive search'),
'example_value': 'true'
}
})
def features(self, request, type_names, property_name, cql_filter=None,
filter_property_name=None, q=None, **kwargs):
params = {
'VERSION': self.get_wfs_service_version(),
'SERVICE': 'WFS',
'REQUEST': 'GetFeature',
'TYPENAMES': type_names,
'PROPERTYNAME': property_name,
'OUTPUTFORMAT': 'json',
}
if cql_filter:
params.update({'CQL_FILTER': cql_filter})
if filter_property_name and q:
if 'case-insensitive' in kwargs:
operator = 'ILIKE'
else:
operator = 'LIKE'
params['CQL_FILTER'] += ' AND %s %s \'%%%s%%\'' % (filter_property_name, operator, q)
response = self.requests.get(self.wfs_service_url, params=params)
data = []
try:
response.json()
except ValueError:
self.handle_opengis_error(response)
# if handle_opengis_error did not raise an error, we raise a generic one
raise APIError(u'OpenGIS Error: unparsable error',
data={'content': repr(response.content[:1024])})
for feature in response.json()['features']:
feature['text'] = feature['properties'].get(property_name)
data.append(feature)
return {'data': data}
def handle_opengis_error(self, response):
try:
root = ET.fromstring(response.content)
except ET.ParseError:
return None
if root.tag != '{http://www.opengis.net/ows/1.1}ExceptionReport':
return None
exception = root.find('{http://www.opengis.net/ows/1.1}Exception')
exception_code = exception.attrib.get('exceptionCode')
if exception is None:
return None
exception_text = exception.find('{http://www.opengis.net/ows/1.1}ExceptionText')
if exception_text is None:
return None
content = exception_text.text
htmlparser = HTMLParser()
content = htmlparser.unescape(content)
raise APIError(u'OpenGIS Error: %s' % exception_code or 'unknown code',
data={'text': content})
def convert_coordinates(self, lon, lat, reverse=False):
lon, lat = float(lon), float(lat)
if self.projection != 'EPSG:4326':
wgs84 = pyproj.Proj(init='EPSG:4326')
target_projection = pyproj.Proj(init=self.projection)
if reverse:
lon, lat = pyproj.transform(target_projection, wgs84, lon, lat)
else:
lon, lat = pyproj.transform(wgs84, target_projection, lon, lat)
return lon, lat
def get_bbox(self, lon1, lat1, lon2, lat2):
if self.projection == 'EPSG:4326':
# send as is but invert coordinates
return '%s,%s,%s,%s' % (lat1, lon1, lat2, lon2)
wgs84 = pyproj.Proj(init='EPSG:4326')
target_projection = pyproj.Proj(init=self.projection)
x1, y1 = pyproj.transform(wgs84, target_projection, lon1, lat1)
x2, y2 = pyproj.transform(wgs84, target_projection, lon2, lat2)
return '%s,%s,%s,%s' % (x1, y1, x2, y2)
@endpoint(perm='can_access',
description=_('Get feature info'),
parameters={
'lat': {'description': _('Latitude'), 'example_value': '45.79689'},
'lon': {'description': _('Longitude'), 'example_value': '4.78414'},
})
def feature_info(self, request, lat, lon):
lat, lon = float(lat), float(lon)
bbox = self.get_bbox(lon - 0.001, lat - 0.001, lon + 0.001, lat + 0.001)
params = {
'VERSION': '1.3.0',
'SERVICE': 'WMS',
'REQUEST': 'GetFeatureInfo',
'INFO_FORMAT': 'application/vnd.ogc.gml',
'STYLES': '',
'I': '0',
'J': '0',
'HEIGHT': '20',
'WIDTH': '20',
'CRS': self.projection,
'LAYERS': self.query_layer,
'QUERY_LAYERS': self.query_layer,
'BBOX': bbox,
}
response = self.requests.get(self.wms_service_url, params=params)
element = ET.fromstring(response.content)
return {'err': 0, 'data': build_dict_from_xml(element)}
# https://carton.entrouvert.org/hydda-tiles/16/33650/23378.pn
@endpoint(name='tile', pattern=r'^(?P<zoom>\d+)/(?P<tile_x>\d+)/(?P<tile_y>\d+).png',
example_pattern='{zoom}/{tile_x}/{tile_y}.png',
parameters={
'zoom': {'description': _('Zoom Level'), 'example_value': '16'},
'tile_x': {'description': _('X Coordinate'), 'example_value': '33650'},
'tile_y': {'description': _('Y Coordinate'), 'example_value': '23378'},
})
def tile(self, request, zoom, tile_x, tile_y):
zoom = int(zoom)
tile_x = int(tile_x)
tile_y = int(tile_y)
def num2deg(xtile, ytile, zoom):
# http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python
n = 2.0 ** zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
lat_deg = math.degrees(lat_rad)
return (lon_deg, lat_deg)
# lower left
ll_lon, ll_lat = num2deg(tile_x, tile_y+1, zoom)
# upper right
ur_lon, ur_lat = num2deg(tile_x+1, tile_y, zoom)
bbox = self.get_bbox(ll_lon, ll_lat, ur_lon, ur_lat)
params = {
'VERSION': '1.3.0',
'SERVICE': 'WMS',
'REQUEST': 'GetMap',
'LAYERS': self.query_layer,
'STYLES': '',
'FORMAT': 'image/png',
'TRANSPARENT': 'false',
'HEIGHT': '256',
'WIDTH': '256',
'CRS': self.projection,
'BBOX': bbox,
}
response = self.requests.get(
self.wms_service_url,
params=params,
cache_duration=300)
return HttpResponse(response.content, content_type='image/png')
@endpoint(perm='can_access', description=_('Get feature info'))
def reverse(self, request, lat, lon, **kwargs):
lon, lat = self.convert_coordinates(lon, lat)
cql_filter = 'DWITHIN(the_geom,Point(%s %s),%s,meters)' % (lon, lat, self.search_radius)
params = {
'VERSION': self.get_wfs_service_version(),
'SERVICE': 'WFS',
'REQUEST': 'GetFeature',
'TYPENAMES': self.query_layer,
'OUTPUTFORMAT': 'json',
'CQL_FILTER': cql_filter
}
response = self.requests.get(self.wfs_service_url, params=params)
if not response.ok:
raise APIError('Webservice returned status code %s' % response.status_code)
closest_feature = {}
min_delta = None
for feature in response.json().get('features'):
if not feature['geometry']['type'] == 'Point':
continue # skip unknown
lon_diff = abs(float(lon) - float(feature['geometry']['coordinates'][0]))
lat_diff = abs(float(lat) - float(feature['geometry']['coordinates'][1]))
# compute hypotenuse til the point
delta = math.sqrt(lon_diff * lon_diff + lat_diff * lat_diff)
if min_delta is None:
min_delta = delta
# choose the shortest
if delta <= min_delta:
closest_feature = feature
if closest_feature:
result = {}
point_lon = closest_feature['geometry']['coordinates'][0]
point_lat = closest_feature['geometry']['coordinates'][1]
point_lon, point_lat = self.convert_coordinates(point_lon, point_lat, reverse=True)
result['lon'] = str(point_lon)
result['lat'] = str(point_lat)
result['address'] = {}
for attribute, properties in self.attributes_mapping:
for field in properties:
if closest_feature['properties'].get(field):
result['address'][attribute] = six.text_type(closest_feature['properties'][field])
break
return result
raise APIError('Unable to geocode')