simplify numpy call

Author: adelepoulle

src/py_eddy_tracker/grid/__init__.py

@@ -3,7 +3,9 @@
 from scipy import interpolate
 from scipy import spatial
 from pyproj import Proj
-import numpy as np
+from numpy import unique, array, unravel_index, r_, floor, interp, arange, \
+    sin, cos, deg2rad, arctan2, sqrt, pi, zeros, reciprocal, ma, empty, \
+    concatenate
 import logging
 from ..tracking_objects import nearest
 
@@ -106,8 +108,8 @@ def read_nc_att(self, varname, att):
 
     @property
     def is_regular(self):
-        steps_lon = np.unique(self._lon[0, 1:] - self._lon[0,:-1])
-        steps_lat = np.unique(self._lat[1:, 0] - self._lat[:-1, 0])
+        steps_lon = unique(self._lon[0, 1:] - self._lon[0,:-1])
+        steps_lat = unique(self._lat[1:, 0] - self._lat[:-1, 0])
         return len(steps_lon) == 1 and len(steps_lat) == 1 and \
             steps_lon[0] != 0. and steps_lat[0] != 0.
 
@@ -135,14 +137,14 @@ def kdt(lon, lat, limits, k=4):
             
             Don't use cKDTree for regular grid
             """
-            ppoints = np.array([lon.ravel(), lat.ravel()]).T
+            ppoints = array([lon.ravel(), lat.ravel()]).T
             ptree = spatial.cKDTree(ppoints)
             pindices = ptree.query(limits, k=k)[1]
-            iind, jind = np.array([], dtype=int), np.array([], dtype=int)
+            iind, jind = array([], dtype=int), array([], dtype=int)
             for pind in pindices.ravel():
-                j, i = np.unravel_index(pind, lon.shape)
-                iind = np.r_[iind, i]
-                jind = np.r_[jind, j]
+                j, i = unravel_index(pind, lon.shape)
+                iind = r_[iind, i]
+                jind = r_[jind, j]
             return iind, jind
 
         if 'AvisoGrid' in self.__class__.__name__:
@@ -152,16 +154,16 @@ def kdt(lon, lat, limits, k=4):
                 Used for a zero crossing, e.g., across Agulhas region
                 """
                 if self.is_regular:
-                    i_1 = int(np.floor(np.interp((lonmin - 0.5) % 360,
+                    i_1 = int(floor(interp((lonmin - 0.5) % 360,
                                        self._lon[0],
-                                       np.arange(len(self._lon[0])))))
-                    i_0 = int(np.floor(np.interp((lonmax + 0.5) % 360,
+                                       arange(len(self._lon[0])))))
+                    i_0 = int(floor(interp((lonmax + 0.5) % 360,
                                        self._lon[0],
-                                       np.arange(len(self._lon[0])))
+                                       arange(len(self._lon[0])))
                                        ) + 1)
                 else:
                     def half_limits(lon, lat):
-                        return np.array([[lon.min(), lon.max(),
+                        return array([[lon.min(), lon.max(),
                                           lon.max(), lon.min()],
                                          [lat.min(), lat.min(),
                                           lat.max(), lat.max()]]).T
@@ -230,12 +232,12 @@ def haversine_dist(self, lon1, lat1, lon2, lat2):
         Return:
           distance (m)
         """
-        sin_dlat = np.sin(np.deg2rad(lat2 - lat1) * 0.5)
-        sin_dlon = np.sin(np.deg2rad(lon2 - lon1) * 0.5)
-        cos_lat1 = np.cos(np.deg2rad(lat1))
-        cos_lat2 = np.cos(np.deg2rad(lat2))
+        sin_dlat = sin(deg2rad(lat2 - lat1) * 0.5)
+        sin_dlon = sin(deg2rad(lon2 - lon1) * 0.5)
+        cos_lat1 = cos(deg2rad(lat1))
+        cos_lat2 = cos(deg2rad(lat2))
         a_val = sin_dlon ** 2 * cos_lat1 * cos_lat2 + sin_dlat ** 2
-        c_val = 2 * np.arctan2(np.sqrt(a_val), np.sqrt(1 - a_val))
+        c_val = 2 * arctan2(sqrt(a_val), sqrt(1 - a_val))
         return 6371315.0 * c_val  # Return the distance
 
     def nearest_point(self, lon, lat):
@@ -260,7 +262,7 @@ def get_aviso_f_pm_pn(self):
         logging.info('--- Computing Coriolis (f), d_x(p_m),'
                      'd_y (p_n) for padded grid')
         # Get GRAVITY / Coriolis
-        self._gof = np.sin(np.deg2rad(self.latpad)) * 4. * np.pi / 86400.
+        self._gof = sin(deg2rad(self.latpad)) * 4. * pi / 86400.
         self._f_val = self._gof.copy()
         self._gof = self.GRAVITY / self._gof
 
@@ -270,29 +272,29 @@ def get_aviso_f_pm_pn(self):
         latv = self.half_interp(self.latpad[:-1], self.latpad[1:])
 
         # Get p_m and p_n
-        p_m = np.zeros_like(self.lonpad)
+        p_m = zeros(self.lonpad.shape)
         p_m[:, 1:-1] = self.haversine_dist(lonu[:, :-1], latu[:, :-1],
                                            lonu[:, 1:], latu[:, 1:])
         p_m[:, 0] = p_m[:, 1]
         p_m[:, -1] = p_m[:, -2]
         self._dx = p_m
-        self._pm = np.reciprocal(p_m)
+        self._pm = reciprocal(p_m)
 
-        p_n = np.zeros_like(self.lonpad)
+        p_n = zeros(self.lonpad.shape)
         p_n[1:-1] = self.haversine_dist(lonv[:-1], latv[:-1],
                                         lonv[1:], latv[1:])
         p_n[0] = p_n[1]
         p_n[-1] = p_n[-2]
         self._dy = p_n
-        self._pn = np.reciprocal(p_n)
+        self._pn = reciprocal(p_n)
         return self
 
     def u2rho_2d(self, uu_in):
         """
         Convert a 2D field at u_val points to a field at rho points
         """
         def uu2ur(uu_in, m_p, l_p):
-            u_out = np.zeros((m_p, l_p))
+            u_out = zeros((m_p, l_p))
             u_out[:, 1:-1] = self.half_interp(uu_in[:, :-1], uu_in[:, 1:])
             u_out[:, 0] = u_out[:, 1]
             u_out[:, -1] = u_out[:, -2]
@@ -303,7 +305,7 @@ def uu2ur(uu_in, m_p, l_p):
     def v2rho_2d(self, vv_in):
         # Convert a 2D field at v_val points to a field at rho points
         def vv2vr(vv_in, m_p, l_p):
-            v_out = np.zeros((m_p, l_p))
+            v_out = zeros((m_p, l_p))
             v_out[1:-1] = self.half_interp(vv_in[:-1], vv_in[1:])
             v_out[0] = v_out[1]
             v_out[-1] = v_out[-2]
@@ -372,13 +374,13 @@ def set_geostrophic_velocity(self, zeta):
         zeta1, zeta2 = zeta.data[1:].view(), zeta.data[:-1].view()
         pn1, pn2 = self.p_n[1:].view(), self.p_n[:-1].view()
         self.upad[:] = self.v2rho_2d(
-            np.ma.array((zeta1 - zeta2) * 0.5 * (pn1 + pn2), mask= self.vmask))
+            ma.array((zeta1 - zeta2) * 0.5 * (pn1 + pn2), mask= self.vmask))
         self.upad *= -self.gof
 
         zeta1, zeta2 = zeta.data[:, 1:].view(), zeta.data[:, :-1].view()
         pm1, pm2 = self.p_m[:, 1:].view(), self.p_m[:, :-1].view()
         self.vpad[:] = self.u2rho_2d(
-            np.ma.array((zeta1 - zeta2) *0.5 * (pm1 + pm2), mask=self.umask))
+            ma.array((zeta1 - zeta2) *0.5 * (pm1 + pm2), mask=self.umask))
         self.vpad *= self.gof
         return self
 
@@ -388,13 +390,13 @@ def set_u_v_eke(self, pad=2):
         """
         j_size = self.slice_j_pad.stop - self.slice_j_pad.start
         if self.zero_crossing:
-            u_1 = np.empty((j_size, self.slice_i_pad.start))
-            u_0 = np.empty((j_size, self._lon.shape[1] - self.slice_i_pad.stop))
-            self.upad = np.ma.concatenate((u_0, u_1), axis=1)
+            u_1 = empty((j_size, self.slice_i_pad.start))
+            u_0 = empty((j_size, self._lon.shape[1] - self.slice_i_pad.stop))
+            self.upad = ma.concatenate((u_0, u_1), axis=1)
         else:
-            self.upad = np.empty((j_size,
+            self.upad = empty((j_size,
                                   self.slice_i_pad.stop - self.slice_i_pad.start))
-        self.vpad = np.empty_like(self.upad)
+        self.vpad = empty(self.upad.shape)
 
     def get_eke(self):
         """
@@ -412,7 +414,7 @@ def uspd(self):
         if hasattr(uspd, 'mask'):
             uspd.mask += self.mask[self.view_unpad]
         else:
-            uspd = np.ma.array(uspd, mask=self.mask[self.view_unpad])
+            uspd = ma.array(uspd, mask=self.mask[self.view_unpad])
         return uspd
 
     def set_interp_coeffs(self, sla, uspd):
@@ -429,8 +431,8 @@ def set_interp_coeffs(self, sla, uspd):
     def create_index_inverse(slice_to_inverse, size):
         """Return an array of index
         """
-        index = np.concatenate((np.arange(slice_to_inverse.stop, size),
-                                np.arange(slice_to_inverse.start)))
+        index = concatenate((arange(slice_to_inverse.stop, size),
+                                arange(slice_to_inverse.start)))
         return index
 
     def gaussian_resolution(self, zwl, mwl):

src/py_eddy_tracker/grid/aviso.py

@@ -3,7 +3,8 @@
 from scipy import ndimage
 from scipy import spatial
 from dateutil import parser
-import numpy as np
+from numpy import meshgrid, zeros, array, where, ma, argmin, vstack, ones \
+    newaxis, sqrt, diff, r_
 import logging
 from netCDF4 import Dataset
 
@@ -15,6 +16,11 @@ class AvisoGrid(PyEddyTracker):
     Class to satisfy the need of the eddy tracker
     to have a grid class
     """
+    KNOWN_UNITS = dict(
+            m=100.,
+            cm=1.,
+            )
+
     def __init__(self, aviso_file, the_domain,
                  lonmin, lonmax, latmin, latmax, grid_name, lon_name,
                  lat_name,with_pad=True):
@@ -41,8 +47,8 @@ def __init__(self, aviso_file, the_domain,
 
         if lonmin < 0 and lonmax <= 0:
             self._lon -= 360.
-        self._lon, self._lat = np.meshgrid(self._lon, self._lat)
-        self._angle = np.zeros_like(self._lon)
+        self._lon, self._lat = meshgrid(self._lon, self._lat)
+        self._angle = zeros(self._lon.shape)
 
         if 'MedSea' in self.the_domain:
             self._lon -= 360.
@@ -71,23 +77,19 @@ def get_aviso_data(self, aviso_file):
         """
         Read nc data from AVISO file
         """
-        KNOWN_UNITS = dict(
-            m=100.,
-            cm=1.,
-            )
         self.grid_filename = aviso_file
         units = self.read_nc_att(self.grid_name, 'units')
-        if units not in KNOWN_UNITS:
+        if units not in self.KNOWN_UNITS:
             raise Exception('Unknown units : %s' % units)
 
         with Dataset(self.grid_filename) as h_nc:
-            grid_dims = np.array(h_nc.variables[self.grid_name].dimensions)
+            grid_dims = array(h_nc.variables[self.grid_name].dimensions)
             lat_dim = h_nc.variables[self.lat_name].dimensions[0]
             lon_dim = h_nc.variables[self.lon_name].dimensions[0]
 
         i_list = []
         transpose = False
-        if np.where(grid_dims == lat_dim)[0][0] > np.where(grid_dims == lon_dim)[0][0]:
+        if where(grid_dims == lat_dim)[0][0] > where(grid_dims == lon_dim)[0][0]:
             transpose = True
         for grid_dim in grid_dims:
             if grid_dim == lat_dim:
@@ -101,11 +103,11 @@ def get_aviso_data(self, aviso_file):
         if transpose:
             zeta = zeta.T
 
-        zeta *= KNOWN_UNITS[units]  # units to cm
+        zeta *= self.KNOWN_UNITS[units]  # units to cm
         if hasattr(zeta, 'mask'):
             return zeta
         else:
-            return np.ma.array(zeta)
+            return ma.array(zeta)
 
     def set_mask(self, sla):
         """
@@ -117,7 +119,7 @@ def set_mask(self, sla):
             if 'Global' in self.the_domain:
 
                 # Close Drake Passage
-                minus70 = np.argmin(np.abs(self.lonpad[0] + 70))
+                minus70 = argmin(abs(self.lonpad[0] + 70))
                 self.mask[:125, minus70] = True
 
                 # DT10 mask is open around Panama, so close it...
@@ -150,8 +152,8 @@ def set_mask(self, sla):
                 self.labels = ndimage.label(-self.mask)[0]
 
                 # Set to known sea point
-                plus200 = np.argmin(np.abs(self.lonpad[0] - 200))
-                plus9 = np.argmin(np.abs(self.latpad[:, 0] - 9))
+                plus200 = argmin(abs(self.lonpad[0] - 200))
+                plus9 = argmin(abs(self.latpad[:, 0] - 9))
                 sea_label = self.labels[plus9, plus200]
                 self.mask += self.labels != sea_label
         return self
@@ -168,8 +170,8 @@ def fillmask(self, data, mask):
 
         # Create (i, j) point arrays for good and bad data.
         # Bad data are marked by the fill_value, good data elsewhere.
-        igood = np.vstack(np.where(data != fill_value)).T
-        ibad = np.vstack(np.where(data == fill_value)).T
+        igood = vstack(where(data != fill_value)).T
+        ibad = vstack(where(data == fill_value)).T
 
         # Create a tree for the bad points, the points to be filled
         tree = spatial.cKDTree(igood)
@@ -180,14 +182,14 @@ def fillmask(self, data, mask):
 
         # Create a normalised weight, the nearest points are weighted as 1.
         #   Points greater than one are then set to zero
-        weight = dist / (dist.min(axis=1)[:, np.newaxis])
-        weight *= np.ones_like(dist)
-        np.place(weight, weight > 1., 0.)
+        weight = dist / (dist.min(axis=1)[:, newaxis])
+        weight *= ones(dist.shape)
+        weight[weight > 1.] = 0.
 
         # Multiply the queried good points by the weight, selecting only the
         # nearest points. Divide by the number of nearest points to get average
         xfill = weight * data[igood[:, 0][iquery], igood[:, 1][iquery]]
-        xfill = (xfill / weight.sum(axis=1)[:, np.newaxis]).sum(axis=1)
+        xfill = (xfill / weight.sum(axis=1)[:, newaxis]).sum(axis=1)
 
         # Place average of nearest good points, xfill, into bad point locations
         data[ibad[:, 0], ibad[:, 1]] = xfill
@@ -262,8 +264,8 @@ def p_n(self):  # Reciprocal of d_y
 
     @property
     def resolution(self):
-        return np.sqrt(np.diff(self.lon[1:], axis=1) *
-                       np.diff(self.lat[:, 1:], axis=0)).mean()
+        return sqrt(diff(self.lon[1:], axis=1) *
+                       diff(self.lat[:, 1:], axis=0)).mean()
 
     @property
     def boundary(self):
@@ -274,8 +276,8 @@ def boundary(self):
         Returns:
           lon/lat boundary points
         """
-        lon = np.r_[(self.lon[:, 0], self.lon[-1],
+        lon = r_[(self.lon[:, 0], self.lon[-1],
                      self.lon[::-1, -1], self.lon[0, ::-1])]
-        lat = np.r_[(self.lat[:, 0], self.lat[-1],
+        lat = r_[(self.lat[:, 0], self.lat[-1],
                      self.lat[::-1, -1], self.lat[0, ::-1])]
         return lon, lat

src/py_eddy_tracker/grid/roms.py

@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-'''
+"""
 ===========================================================================
 This file is part of py-eddy-tracker.
 
@@ -30,9 +30,9 @@
 Useful when the output files don't contain the grid
 information
 
-'''
+"""
+
 import netCDF4 as netcdf
-# from matplotlib.mlab import load
 import numpy as np
 import matplotlib.path as Path
 from . import PyEddyTracker

src/py_eddy_tracker/observations.py

@@ -141,7 +141,7 @@ def merge(self, other):
         return eddies
 
     def reset(self):
-        self.observations = np.zeros(0, dtype=self.dtype)
+        self.observations = zeros(0, dtype=self.dtype)
 
     @property
     def obs(self):