1# (C) British Crown Copyright 2014 - 2018, Met Office
2#
3# This file is part of pyepsg.
4#
5# pyepsg is free software: you can redistribute it and/or modify it under
6# the terms of the GNU Lesser General Public License as published by the
7# Free Software Foundation, either version 3 of the License, or
8# (at your option) any later version.
9#
10# pyepsg is distributed in the hope that it will be useful,
11# but WITHOUT ANY WARRANTY; without even the implied warranty of
12# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13# GNU Lesser General Public License for more details.
14#
15# You should have received a copy of the GNU Lesser General Public License
16# along with pyepsg.  If not, see <http://www.gnu.org/licenses/>.
17"""
18Provides simple access to http://epsg.io/.
19
20The entry point for this package is the :func:`get()` function.
21
22"""
23from __future__ import print_function
24
25import sys
26import weakref
27import xml.etree.ElementTree as ET
28
29import requests
30
31
32EPSG_IO_URL = 'http://epsg.io/'
33
34GML_NS = '{http://www.opengis.net/gml/3.2}'
35XLINK_NS = '{http://www.w3.org/1999/xlink}'
36
37
38_cache = weakref.WeakValueDictionary()
39
40
41class EPSG(object):
42    """Parent class of all objects returned by pyepsg."""
43    def __init__(self, element):
44        self.element = element
45
46    @property
47    def identifier(self):
48        """The official URN for this object."""
49        return self.element.find(GML_NS + 'identifier').text
50
51
52class UOM(EPSG):
53    """A unit of measure."""
54    @property
55    def name(self):
56        """The human-readable name."""
57        return self.element.find(GML_NS + 'name').text
58
59
60class Axis(EPSG):
61    """A single coordinate axis."""
62
63    @property
64    def direction(self):
65        """A description of the orientation of this axis."""
66        return self.element.find(GML_NS + 'axisDirection').text
67
68    @property
69    def uom(self):
70        """The name of the unit of measure used on this axis."""
71        uom = self.element.attrib['uom']
72        return get(uom).name
73
74    def __repr__(self):
75        return '<Axis: {self.direction} / {self.uom}>'.format(self=self)
76
77
78class CartesianCS(EPSG):
79    """A 1-, 2-, or 3-dimensional cartesian coordinate system."""
80    @property
81    def axes(self):
82        """An ordered list of :class:`Axis` objects describing X and Y."""
83        axes = self.element.findall(GML_NS + 'axis')
84        return [Axis(axis.find(GML_NS + 'CoordinateSystemAxis')) for
85                axis in axes]
86
87    @property
88    def name(self):
89        """The human-readable name."""
90        return self.element.find(GML_NS + 'name').text
91
92    @property
93    def remarks(self):
94        """Human-readable comments."""
95        return self.element.find(GML_NS + 'remarks').text
96
97    def __repr__(self):
98        name = self.name
99        if len(name) > 38:
100            name = name[:38] + '..'
101        return '<CartesianCS: {name}>'.format(name=name)
102
103
104class CRS(EPSG):
105    """
106    Abstract parent class for :class:`GeodeticCRS`, :class:`ProjectedCRS`
107    and :class:`CompoundCRS`.
108
109    """
110
111    @property
112    def id(self):
113        """The EPSG code for this CRS."""
114        id = self.element.attrib[GML_NS + 'id']
115        code = id.split('-')[-1]
116        return code
117
118    @property
119    def name(self):
120        """The human-readable name."""
121        return self.element.find(GML_NS + 'name').text
122
123    @property
124    def scope(self):
125        """A human-readable description of the intended usage for this CRS."""
126        return self.element.find(GML_NS + 'scope').text
127
128    def as_esri_wkt(self):
129        """
130        Return the ESRI WKT string which corresponds to the CRS.
131
132        For example::
133
134            >>> print(get(27700).as_esri_wkt())  # doctest: +ELLIPSIS
135            PROJCS["OSGB_1936_British_National_Grid",GEOGCS["GCS_OSGB 19...
136
137        """
138        url = '{prefix}{code}.esriwkt?download'.format(prefix=EPSG_IO_URL,
139                                                       code=self.id)
140        return requests.get(url).text
141
142    def as_html(self):
143        """
144        Return the OGC WKT which corresponds to the CRS as HTML.
145
146        For example::
147
148            >>> print(get(27700).as_html())  # doctest: +ELLIPSIS
149            <div class="syntax"><pre><span class="gh">PROJCS</span><span...
150
151        """
152        url = '{prefix}{code}.html?download'.format(prefix=EPSG_IO_URL,
153                                                    code=self.id)
154        return requests.get(url).text
155
156    def as_proj4(self):
157        """
158        Return the PROJ.4 string which corresponds to the CRS.
159
160        For example::
161
162            >>> print(get(21781).as_proj4())
163            +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 \
164+k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel \
165+towgs84=674.4,15.1,405.3,0,0,0,0 +units=m +no_defs
166
167        """
168        url = '{prefix}{code}.proj4?download'.format(prefix=EPSG_IO_URL,
169                                                     code=self.id)
170        return requests.get(url).text.strip()
171
172    def as_wkt(self):
173        """
174        Return the OGC WKT string which corresponds to the CRS.
175
176        For example::
177
178            >>> print(get(27700).as_wkt())  # doctest: +ELLIPSIS
179            PROJCS["OSGB 1936 / British National Grid",GEOGCS["OSGB 1936...
180
181        """
182        url = '{prefix}{code}.wkt?download'.format(prefix=EPSG_IO_URL,
183                                                   code=self.id)
184        return requests.get(url).text
185
186    def domain_of_validity(self):
187        """
188        Return the domain of validity for this CRS as:
189        (west, east, south, north).
190
191        For example::
192
193            >>> print(get(21781).domain_of_validity())
194            [5.96, 10.49, 45.82, 47.81]
195
196
197        """
198        # TODO: Generalise interface to return a polygon? (Can we find
199        # something that uses a polygon instead?)
200        domain = self.element.find(GML_NS + 'domainOfValidity')
201        domain_href = domain.attrib[XLINK_NS + 'href']
202        url = '{prefix}{code}.gml?download'.format(prefix=EPSG_IO_URL,
203                                                   code=domain_href)
204        xml = requests.get(url).content
205        gml = ET.fromstring(xml)
206
207        def extract_bound(tag):
208            ns = '{http://www.isotc211.org/2005/gmd}'
209            xpath = './/{ns}EX_GeographicBoundingBox/{ns}{tag}/'.format(
210                ns=ns,
211                tag=tag)
212            bound = gml.find(xpath)
213            return float(bound.text)
214
215        tags = ('westBoundLongitude', 'eastBoundLongitude',
216                'southBoundLatitude', 'northBoundLatitude')
217        bounds = [extract_bound(tag) for tag in tags]
218        return bounds
219
220
221class GeodeticCRS(CRS):
222    """
223    Represents a single geodetic CRS.
224
225    """
226
227    def __repr__(self):
228        return '<GeodeticCRS: {self.id}, {self.name}>'.format(self=self)
229
230
231class ProjectedCRS(CRS):
232    """
233    Represents a single projected CRS.
234
235    """
236    @property
237    def base_geodetic_crs(self):
238        """The :class:`GeodeticCRS` on which this projection is based."""
239        base = self.element.find(GML_NS + 'baseGeodeticCRS')
240        href = base.attrib[XLINK_NS + 'href']
241        return get(href)
242
243    @property
244    def cartesian_cs(self):
245        """The :class:`CartesianCS` which describes the coordinate axes."""
246        cs = self.element.find(GML_NS + 'cartesianCS')
247        href = cs.attrib[XLINK_NS + 'href']
248        return get(href)
249
250    def __repr__(self):
251        return '<ProjectedCRS: {self.id}, {self.name}>'.format(self=self)
252
253
254class CompoundCRS(CRS):
255    """
256    Represents a single compound CRS.
257
258    """
259    def __repr__(self):
260        return '<CompoundCRS: {self.id}, {self.name}>'.format(self=self)
261
262
263def get(code):
264    """
265    Return an object that corresponds to the given EPSG code.
266
267    Currently supported object types are:
268        - :class:`GeodeticCRS`
269        - :class:`ProjectedCRS`
270        - :class:`CartesianCS`
271        - :class:`UOM`
272
273    For example::
274
275        >>> print(get(27700))
276        <ProjectedCRS: 27700, OSGB 1936 / British National Grid>
277        >>> print(get('4400-cs'))
278        <CartesianCS: Cartesian 2D CS. Axes: easting, northi..>
279        >>> print(get(5973))
280        <CompoundCRS: 5973, ETRS89 / UTM zone 33 + NN2000 height>
281
282    """
283    instance = _cache.get(code)
284    if instance is None:
285        url = '{prefix}{code}.gml?download'.format(prefix=EPSG_IO_URL,
286                                                   code=code)
287        xml = requests.get(url).content
288        root = ET.fromstring(xml)
289        class_for_tag = {
290            GML_NS + 'CartesianCS': CartesianCS,
291            GML_NS + 'GeodeticCRS': GeodeticCRS,
292            GML_NS + 'ProjectedCRS': ProjectedCRS,
293            GML_NS + 'CompoundCRS': CompoundCRS,
294            GML_NS + 'BaseUnit': UOM,
295        }
296        if root.tag in class_for_tag:
297            instance = class_for_tag[root.tag](root)
298        else:
299            raise ValueError('Unsupported code type: {}'.format(root.tag))
300        _cache[code] = instance
301    return instance
302
303
304if __name__ == '__main__':
305    import doctest
306    failure_count, test_count = doctest.testmod()
307    sys.exit(failure_count)
308