--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/pyams_gis/point.py	Thu May 18 17:23:48 2017 +0200
@@ -0,0 +1,64 @@
+#
+# Copyright (c) 2008-2017 Thierry Florac <tflorac AT ulthar.net>
+# All Rights Reserved.
+#
+# This software is subject to the provisions of the Zope Public License,
+# Version 2.1 (ZPL).  A copy of the ZPL should accompany this distribution.
+# THIS SOFTWARE IS PROVIDED "AS IS" AND ANY AND ALL EXPRESS OR IMPLIED
+# WARRANTIES ARE DISCLAIMED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF TITLE, MERCHANTABILITY, AGAINST INFRINGEMENT, AND FITNESS
+# FOR A PARTICULAR PURPOSE.
+#
+
+__docformat__ = 'restructuredtext'
+
+
+# import standard library
+try:
+    from osgeo.osr import SpatialReference, CoordinateTransformation
+    have_gdal = True
+except ImportError:
+    have_gdal = False
+from persistent import Persistent
+
+# import interfaces
+from pyams_gis.interfaces import IGeoPoint, IGeoPointZ, WGS84
+
+# import packages
+from zope.interface import implementer
+from zope.schema.fieldproperty import FieldProperty
+
+
+@implementer(IGeoPoint)
+class GeoPoint(Persistent):
+    """GeoPoint attribute object"""
+
+    longitude = FieldProperty(IGeoPoint['longitude'])
+    latitude = FieldProperty(IGeoPoint['latitude'])
+    projection = FieldProperty(IGeoPoint['projection'])
+
+    def __bool__(self):
+        return bool(self.longitude and self.latitude)
+
+    def get_coordinates(self, projection=WGS84):
+        if projection == self.projection:
+            return self.longitude, self.latitude
+        if (not have_gdal) or not self:
+            return None, None
+        source = SpatialReference()
+        source.ImportFromEPSG(self.projection)
+        destination = SpatialReference()
+        destination.ImportFromEPSG(projection)
+        transformation = CoordinateTransformation(source, destination)
+        return transformation.TransformPoint(float(self.longitude), float(self.latitude))[0:2]
+
+    @property
+    def wgs_coordinates(self):
+        return self.get_coordinates(WGS84)
+
+
+@implementer(IGeoPointZ)
+class GeoPointZ(GeoPoint):
+    """GeoPointZ attribute object"""
+
+    altitude = FieldProperty(IGeoPointZ['altitude'])