passerelle/passerelle/apps/opengis/models.py

360 lines
15 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
import six
import pyproj
from django.core.cache import cache
from django.db import models
from django.http import HttpResponse
from django.utils.six.moves.html_parser import HTMLParser
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.conversion import num2deg
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_capabilities(self, service_type, service_url):
if not service_url:
raise APIError('no %s URL declared' % service_type)
return self.requests.get(
service_url,
params={'request': 'GetCapabilities', 'service': service_type.upper()},
)
def get_service_version(self, service_type, service_url, renew=False):
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.get_capabilities(service_type, service_url)
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 get_output_format(self):
cache_key = 'opengis-%s-output-format' % self.id
output_format = cache.get(cache_key)
if output_format:
return output_format
response = self.get_capabilities('wfs', self.wfs_service_url)
element = ET.fromstring(response.content)
ns = {'ows': 'http://www.opengis.net/ows/1.1'}
formats = element.findall('.//ows:Operation[@name="GetFeature"]/'
'ows:Parameter[@name="outputFormat"]/'
'ows:AllowedValues/ows:Value', ns)
for output_format in formats:
if 'json' in output_format.text.lower():
cache.set(cache_key, output_format.text, 3600)
return output_format.text
raise APIError('WFS server doesn\'t support json output format for GetFeature request')
def get_typename_label(self):
version_str = self.get_wfs_service_version()
version_tuple = tuple(int(x) for x in version_str.split('.'))
if version_tuple <= (1, 1, 0):
return 'TYPENAME'
else:
return 'TYPENAMES'
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()
def build_get_features_params(self, typename=None, property_name=None, cql_filter=None):
params = {
'version': self.get_wfs_service_version(),
'service': 'WFS',
'request': 'GetFeature',
self.get_typename_label(): typename,
'propertyName': property_name,
'outputFormat': self.get_output_format(),
}
if cql_filter:
params['cql_filter'] = cql_filter
return params
@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):
if cql_filter:
if filter_property_name and q:
if 'case-insensitive' in kwargs:
operator = 'ILIKE'
else:
operator = 'LIKE'
cql_filter += ' AND %s %s \'%%%s%%\'' % (filter_property_name, operator, q)
params = self.build_get_features_params(type_names, property_name, cql_filter)
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 '%.6f,%.6f,%.6f,%.6f' % (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 '%.6f,%.6f,%.6f,%.6f' % (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):
try:
lat, lon = float(lat), float(lon)
except ValueError:
raise APIError('Bad coordinates format')
bbox = self.get_bbox(lon - 0.0001, lat - 0.0001, lon + 0.0001, lat + 0.0001)
params = {
'VERSION': '1.3.0',
'SERVICE': 'WMS',
'REQUEST': 'GetFeatureInfo',
'INFO_FORMAT': 'application/vnd.ogc.gml',
'STYLES': '',
'I': '24', 'J': '24', # pixel in the middle of
'HEIGHT': '50', 'WIDTH': '50', # a 50x50 square
'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)
# 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(%.6f %.6f),%s,meters)' % (lon, lat, self.search_radius)
params = self.build_get_features_params(typename=self.query_layer, 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'] = "%.6f" % point_lon
result['lat'] = "%.6f" % 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')