some small changes so now tests run appropriately

tests/test_grid.py

@@ -143,22 +143,21 @@ def test_interpolation():
     # There is a shift between the rho grid and the grid space grid because
     # of the staggered layout. Grid space counts from the u/v grid and
     # therefore is a little different from the rho grid.
-    assert np.allclose(X, grid.x_rho[2, 3] - 0.5)
-    assert np.allclose(Y, grid.y_rho[2, 3] - 0.5)
+    assert np.allclose(X, grid.x_rho[2, 3] + 0.5)
+    assert np.allclose(Y, grid.y_rho[2, 3] + 0.5)
 
     # grid space to projected coordinates, delaunay
     x, y, _ = tracpy.tools.interpolate2d(grid.X[2, 3], grid.Y[2, 3], grid,
                                          'd_ij2xy')
-
-    assert np.allclose(x, grid.X[2, 3] + 0.5)
-    assert np.allclose(y, grid.Y[2, 3] + 0.5)
+    assert np.allclose(x, grid.X[2, 3] - 0.5)
+    assert np.allclose(y, grid.Y[2, 3] - 0.5)
 
     # grid space to projected coordinates, map_coords
     x, y, _ = tracpy.tools.interpolate2d(grid.X[2, 3], grid.Y[2, 3], grid,
                                          'm_ij2xy')
 
-    assert np.allclose(x, grid.X[2, 3] + 0.5)
-    assert np.allclose(y, grid.Y[2, 3] + 0.5)
+    assert np.allclose(x, grid.X[2, 3] - 0.5)
+    assert np.allclose(y, grid.Y[2, 3] - 0.5)
 
 
 # def test_interpolation():

tests/test_rect.py

@@ -308,3 +308,96 @@ def test_run_2d_xy_ssh():
     print(xp[:, -1] - xp[:, 0])
 
     assert np.allclose(xp[:, -1] - xp[:, 0], distance)
+
+
+def test_run_3d_ll():
+    """
+    Test initialization and running of tracpy with lat/lon coordinates in 3d.
+
+    can we initialize and run tracpy (using rectangle example). Compare final
+    location of drifters with known analytic answer. Using lon/lat coords.
+    """
+
+    # some simple example data
+    currents_filename = os.path.join('input', 'ocean_his_0001.nc')
+    grid_filename = os.path.join('input', 'grid.nc')
+    time_units = 'seconds since 1970-01-01'
+    num_layers = 3
+
+    name = 'test_run_3d_ll'
+
+    # Start date in date time formatting
+    date = datetime.datetime(2013, 12, 19, 0)
+
+    # Time between outputs
+    tseas = 4*3600.  # 4 hours between outputs, in seconds
+
+    # Number of days to run the drifters.
+    ndays = tseas*9./(3600.*24)
+
+    # Sets a smaller limit than between model outputs for when to force
+    # interpolation if hasn't already occurred.
+    nsteps = 5
+
+    # Controls the sampling frequency of the drifter tracks.
+    N = 4
+
+    # This allows the user to call to TRACMASS for a different period of time
+    # than between 2 model outputs
+    # Just testing to try new loop, should have same behavior as before
+    dtFromTracmass = tseas/2.
+
+    # Use ff = 1 for forward in time and ff = -1 for backward in time.
+    ff = 1  # will work for ff=1 or ff=-1 since checks by distance traveled
+
+    ah = 0.  # m^2/s
+    av = 0.  # m^2/s
+
+    # turbulence/diffusion flag
+    doturb = 0
+
+    # two particles (starting positions)
+    lon0 = [-123., -123.]
+    lat0 = [48.55, 48.75]
+
+    do3d = 1  # flag to set to 2-d
+
+    z0 = [0,0]  # 'z' 'salt' 's'
+    zpar = 'fromZeta'  # top layer
+
+    usespherical = True
+
+    # Get projection object
+    proj = tracpy.tools.make_proj(setup='galveston', usebasemap=False)
+
+    # Read in grid
+    grid = tracpy.inout.readgrid(grid_filename, proj,
+                                 vert_filename=currents_filename,
+                                 usespherical=usespherical)
+
+    # Initialize Tracpy class
+    tp = Tracpy(currents_filename, grid, name=name, tseas=tseas,
+                ndays=ndays, nsteps=nsteps, N=N, ff=ff, ah=ah, av=av,
+                doturb=doturb, do3d=do3d, z0=z0, zpar=zpar,
+                time_units=time_units, dtFromTracmass=dtFromTracmass,
+                usespherical=usespherical)
+
+    lonp, latp, zp, t, T0, U, V = tracpy.run.run(tp, date, lon0, lat0)
+
+    ## check the results:
+    print(lonp.shape)
+    print(lonp)
+    print(latp)
+
+    # since eastward current, latitude should not change:
+    assert np.allclose(lat0, latp.T)
+
+    # current velocity -- 0.1 m/s
+    # position
+    distance = (ndays * 24 * 3600 * 0.1)*ff
+
+    # better to use pyproj to compute the geodesic
+    geod = pyproj.Geod(ellps='clrk66')
+    end = geod.fwd(lon0, lat0, (90, 90), (distance, distance), radians=False)
+
+    assert np.allclose(lonp[:, -1], end[0])

tracpy/tracpy_class.py

@@ -229,8 +229,9 @@ def prepare_for_model_run(self, date, lon0, lat0):
         # Do z a little lower down
 
         # Initialize seed locations
-        ia = np.ceil(xstart0)
-        ja = np.ceil(ystart0)
+        # these will be used as indices so must be ints
+        ia = np.ceil(xstart0).astype(int)
+        ja = np.ceil(ystart0).astype(int)
 
         # don't use nan's
         # pdb.set_trace()
@@ -309,7 +310,7 @@ def prepare_for_model_run(self, date, lon0, lat0):
         else:   # 3d case
             # Convert initial real space vertical locations to grid space
             # first find indices of grid cells vertically
-            ka = np.ones(ia.size)*np.nan
+            ka = np.ones(ia.size, dtype=int)*-999  # need int placeholder
             zstart0 = np.ones(ia.size)*np.nan
 
             if self.zpar == 'fromMSL':