Merge f3914b9 into 7b89abe

nupic/encoders/geospatial_coordinate.py

@@ -22,15 +22,17 @@
 import math
 
 import numpy
-from pyproj import Proj
+from pyproj import Proj, transform
 from nupic.encoders.coordinate import CoordinateEncoder
 
 
 
-PROJ = Proj(init="epsg:3785")  # Spherical Mercator
 # From http://spatialreference.org/ref/epsg/popular-visualisation-crs-mercator/
+PROJ = Proj(init="epsg:3785")  # Spherical Mercator
 PROJ_RANGE=(20037508.3428, 19971868.8804)  # in meters
 
+# See http://gis.stackexchange.com/a/73829/41082
+geocentric = Proj('+proj=geocent +datum=WGS84 +units=m +no_defs')
 
 
 class GeospatialCoordinateEncoder(CoordinateEncoder):
@@ -66,34 +68,44 @@ def __init__(self,
 
   def getDescription(self):
     """See `nupic.encoders.base.Encoder` for more information."""
-    return [('longitude', 0), ('latitude', 1), ('speed', 2)]
+    return [('speed', 0), ('longitude', 1), ('latitude', 2), ('altitude', 3)]
 
 
   def encodeIntoArray(self, inputData, output):
     """
     See `nupic.encoders.base.Encoder` for more information.
 
-    @param inputData (tuple) Contains longitude (float),
-                             latitude (float), speed (float)
+    @param inputData (tuple) Contains speed (float), longitude (float),
+                             latitude (float), altitude (float)
     @param output (numpy.array) Stores encoded SDR in this numpy array
     """
-    (longitude, latitude, speed) = inputData
-    coordinate = self.coordinateForPosition(longitude, latitude)
+    altitude = None
+    if len(inputData) == 4:
+      (speed, longitude, latitude, altitude) = inputData
+    else:
+      (speed, longitude, latitude) = inputData
+    coordinate = self.coordinateForPosition(longitude, latitude, altitude)
     radius = self.radiusForSpeed(speed)
     super(GeospatialCoordinateEncoder, self).encodeIntoArray(
      (coordinate, radius), output)
 
 
-  def coordinateForPosition(self, longitude, latitude):
+  def coordinateForPosition(self, longitude, latitude, altitude=None):
     """
     Returns coordinate for given GPS position.
 
     @param longitude (float) Longitude of position
     @param latitude (float) Latitude of position
+    @param altitude (float) Altitude of position
     @return (numpy.array) Coordinate that the given GPS position
                           maps to
     """
-    coordinate = numpy.array(PROJ(longitude, latitude))
+    coords = PROJ(longitude, latitude)
+
+    if altitude is not None:
+      coords = transform(PROJ,geocentric,coords[0],coords[1],altitude)
+
+    coordinate = numpy.array(coords)
     coordinate = coordinate / self.scale
     return coordinate.astype(int)
 

tests/unit/nupic/encoders/geospatial_coordinate_test.py

@@ -44,6 +44,25 @@ def testCoordinateForPosition(self):
     self.assertEqual(coordinate.tolist(), [-453549, 150239])
 
 
+  def testCoordinateForPosition3D(self):
+    scale = 30 # meters
+    encoder = GeospatialCoordinateEncoder(scale, 60)
+    coordinate = encoder.coordinateForPosition(
+      -122.229194, 37.486782, 1500
+    )
+    self.assertEqual(coordinate.tolist(), [-90102, -142918, 128710])
+
+
+  def testCoordinateForPositionOrigin3D(self):
+    scale = 1 # meters
+    encoder = GeospatialCoordinateEncoder(scale, 60)
+    coordinate = encoder.coordinateForPosition(0,0,0)
+
+    # see WGS80 defining parameters (semi-major axis) on
+    # http://en.wikipedia.org/wiki/Geodetic_datum#Parameters_for_some_geodetic_systems
+    self.assertEqual(coordinate.tolist(), [6378137, 0, 0])  
+
+
   def testCoordinateForPositionOrigin(self):
     scale = 30  # meters
     encoder = GeospatialCoordinateEncoder(scale, 60)
@@ -88,20 +107,54 @@ def testEncodeIntoArray(self):
     encoder = GeospatialCoordinateEncoder(scale, timestep,
                                           n=999,
                                           w=25)
-    encoding1 = encode(encoder, -122.229194, 37.486782, speed)
-    encoding2 = encode(encoder, -122.229294, 37.486882, speed)
-    encoding3 = encode(encoder, -122.229294, 37.486982, speed)
+    encoding1 = encode(encoder, speed, -122.229194, 37.486782)
+    encoding2 = encode(encoder, speed, -122.229294, 37.486882)
+    encoding3 = encode(encoder, speed, -122.229294, 37.486982)
+
+    overlap1 = overlap(encoding1, encoding2)
+    overlap2 = overlap(encoding1, encoding3)
+
+    self.assertTrue(overlap1 > overlap2)
+
+
+  def testEncodeIntoArrayAltitude(self):
+    scale = 30 # meters
+    timestep = 60 # seconds
+    speed = 2.5 # meters per second
+    longitude, latitude = -122.229294, 37.486782
+    encoder = GeospatialCoordinateEncoder(scale, timestep,
+                                          n=999,
+                                          w=25)
+    encoding1 = encode(encoder, speed, longitude, latitude, 0)
+    encoding2 = encode(encoder, speed, longitude, latitude, 100)
+    encoding3 = encode(encoder, speed, longitude, latitude, 1000)
 
     overlap1 = overlap(encoding1, encoding2)
     overlap2 = overlap(encoding1, encoding3)
 
     self.assertTrue(overlap1 > overlap2)
 
 
+  def testEncodeIntoArray3D(self):
+    scale = 30 # meters
+    timestep = 60 # seconds
+    speed = 2.5 # meters per second
+    encoder = GeospatialCoordinateEncoder(scale, timestep,
+                                          n=999,
+                                          w=25)
+    encoding1 = encode(encoder, speed, -122.229194, 37.486782, 0)
+    encoding2 = encode(encoder, speed, -122.229294, 37.486882, 100)
+    encoding3 = encode(encoder, speed, -122.229294, 37.486982, 1000)
+
+    overlap1 = overlap(encoding1, encoding2)
+    overlap2 = overlap(encoding1, encoding3)
+
+    self.assertTrue(overlap1 > overlap2)
+
 
-def encode(encoder, longitude, latitude, speed):
+def encode(encoder, speed, longitude, latitude, altitude=None):
   output = np.zeros(encoder.getWidth(), dtype=defaultDtype)
-  encoder.encodeIntoArray((longitude, latitude, speed), output)
+  encoder.encodeIntoArray((speed, longitude, latitude, altitude), output)
   return output
 
 def overlap(sdr1, sdr2):