Repository: ClimateImpactLab/climate_toolbox

Issue:
vectorized distance functions
something along these lines... (some code courtesy of [geopy](https://github.com/geopy/geopy/blob/master/geopy/distance.py)):

```python

#             model             major (km)   minor (km)     flattening
ELLIPSOIDS = {'WGS-84':        (6378.137,    6356.7523142,  1 / 298.257223563),
              'GRS-80':        (6378.137,    6356.7523141,  1 / 298.257222101),
              'Airy (1830)':   (6377.563396, 6356.256909,   1 / 299.3249646),
              'Intl 1924':     (6378.388,    6356.911946,   1 / 297.0),
              'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465),
              'GRS-67':        (6378.1600,   6356.774719,   1 / 298.25),
              }


EARTH_RADIUS = 6371.009

def great_circle(ax, ay, bx, by, radius=EARTH_RADIUS):
    '''
    calculate the great circle distance (km) between points
    
    Provide points (ax, ay) and (bx, by) as floats, or as 
    vectors. If ax and ay are vectors or arrays of the
    same shape, the element-wise distance will be found
    between points in the vectors/arrays. If ax, ay are
    (Mx1) column vectors and (bx, by) are (1xN) row 
    vectors, the vectors will be broadcast using numpy 
    broadcasting rules and the distance between each pair
    of points will be returned as an (MxN) matrix.

    Parameters
    -----------
    ax : float or array
        x/long of point a
    ay : float or array
        y/lat of point a
    bx : float or array
        x/long of point b
    by : float or array
        y/lat of point b
    radius : float, optional
        Radius of the sphere on which to calculate the great
        circle distance (default is to use the Earth's radius in
        km, `6371.009`). Values returned will be in units of the
        radius provided. 

    Returns
    --------
    distance : float or array
        great circle distance between points a and b. Units will
        match the radius provided (default km)
    '''

    lat1, lng1 = np.radians(ay), np.radians(ax)
    lat2, lng2 = np.radians(by), np.radians(bx)

    sin_lat1, cos_lat1 = np.sin(lat1), np.cos(lat1)
    sin_lat2, cos_lat2 = np.sin(lat2), np.cos(lat2)

    delta_lng = lng2 - lng1
    cos_delta_lng, sin_delta_lng = np.cos(delta_lng), np.sin(delta_lng)

    d = np.arctan2(np.sqrt((cos_lat2 * sin_delta_lng) ** 2 +
                   (cos_lat1 * sin_lat2 -
                    sin_lat1 * cos_lat2 * cos_delta_lng) ** 2),
              sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)

    return radius * d
```

# todo

add [geodesic](https://github.com/geopy/geopy/blob/master/geopy/distance.py#L354) distance

Discussion:

Patch:
diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
index 474c8a5..5d223c0 100644
--- a/CONTRIBUTING.rst
+++ b/CONTRIBUTING.rst
@@ -70,13 +70,6 @@ Ready to contribute? Here's how to set up `climate_toolbox` for local developmen
     $ cd climate_toolbox/
     $ python setup.py develop
 
-   If you are using Conda run the following::
-
-    $ conda create -n climate_toolbox python=3.5
-    $ conda activate climate_toolbox
-    $ pip install -r requirements.txt
-    $ python setup.py test 
-
 4. Create a branch for local development::
 
     $ git checkout -b name-of-your-bugfix-or-feature
diff --git a/climate_toolbox/aggregations/aggregations.py b/climate_toolbox/aggregations/aggregations.py
deleted file mode 100644
index bbaf97d..0000000
--- a/climate_toolbox/aggregations/aggregations.py
+++ /dev/null
@@ -1,146 +0,0 @@
-import xarray as xr
-import numpy as np
-import pandas as pd
-
-from distutils.version import LooseVersion
-
-
-def _reindex_spatial_data_to_regions(ds, df):
-    '''
-    Reindexes spatial and segment weight data to regions
-    Enables region index-based math operations
-    Parameters
-    ----------
-    ds: xarray Dataset
-    df: pandas DataFrame
-    Returns
-    -------
-    Xarray DataArray
-    '''
-
-    # use vectorized indexing in xarray >= 0.10
-    if LooseVersion(xr.__version__) > LooseVersion('0.9.999'):
-
-        lon_indexer = xr.DataArray(df.lon.values, dims=('reshape_index', ))
-        lat_indexer = xr.DataArray(df.lat.values, dims=('reshape_index', ))
-
-        return ds.sel(lon=lon_indexer, lat=lat_indexer)
-
-    else:
-        res = ds.sel_points(
-            'reshape_index',
-            lat=df.lat.values,
-            lon=df.lon.values)
-
-        return res
-
-
-def _aggregate_reindexed_data_to_regions(
-        ds,
-        variable,
-        aggwt,
-        agglev,
-        weights,
-        backup_aggwt='areawt'):
-    '''
-    Performs weighted avg for climate variable by region
-
-    Parameters
-    ----------
-
-    ds: xarray.DataArray
-
-    variable: str
-        name of the data variable
-
-    aggwt: str
-        variable to weight by (i.e popwt, areawt, cropwt)
-
-    agglev: str
-        indicates which regional id scheme to select in the dataframe
-
-    weight: pd.DataFrame
-        pandas DataFrame of weights
-
-    backup_aggwt: str, optional
-        aggregation weight to use in regions with no aggwt data (default
-        'areawt')
-
-    '''
-
-    ds.coords[agglev] = xr.DataArray(
-                weights[agglev].values,
-                dims={'reshape_index': weights.index.values})
-
-    # format weights
-    ds[aggwt] = xr.DataArray(
-                weights[aggwt].values,
-                dims={'reshape_index': weights.index.values})
-
-    ds[aggwt] = (
-        ds[aggwt]
-        .where(ds[aggwt] > 0)
-        .fillna(weights[backup_aggwt].values))
-
-    weighted = xr.Dataset({
-        variable: (
-            (
-                (ds[variable]*ds[aggwt])
-                .groupby(agglev)
-                .sum(dim='reshape_index')) /
-            (
-                ds[aggwt]
-                .groupby(agglev)
-                .sum(dim='reshape_index')))})
-
-    return weighted
-
-
-def weighted_aggregate_grid_to_regions(
-        ds,
-        variable,
-        aggwt,
-        agglev,
-        weights=None):
-    '''
-    Computes the weighted reshape of gridded data
-
-    Parameters
-    ----------
-    ds : xr.Dataset
-        xarray Dataset to be aggregated. Must have 'lat' and 'lon' in the
-        coordinates.
-
-    variable : str
-        name of the variable to be aggregated
-
-    aggwt : str
-        Weighting variable (e.g. 'popwt', 'areawt'). This must be a column name
-        in the weights file.
-
-    agglev : str
-        Target regional aggregation level (e.g. 'ISO', 'hierid'). This must be
-        a column name in the weights file.
-
-    weights : str, optional
-        Regional aggregation weights (default agglomerated-world-new BCSD
-        segment weights)
-
-    Returns
-    -------
-    ds: xr.Dataset
-        weighted and averaged dataset based on agglev
-    '''
-
-    if weights is None:
-        weights = _prepare_spatial_weights_data()
-
-    ds = _reindex_spatial_data_to_regions(ds, weights)
-    ds = _aggregate_reindexed_data_to_regions(
-        ds,
-        variable,
-        aggwt,
-        agglev,
-        weights)
-
-    return ds
diff --git a/climate_toolbox/cli.py b/climate_toolbox/cli.py
new file mode 100644
index 0000000..38a125d
--- /dev/null
+++ b/climate_toolbox/cli.py
@@ -0,0 +1,17 @@
+# -*- coding: utf-8 -*-
+
+"""Console script for climate_toolbox."""
+
+import click
+
+
+@click.command()
+def main(args=None):
+    """Console script for climate_toolbox."""
+    click.echo("Replace this message by putting your code into "
+               "climate_toolbox.cli.main")
+    click.echo("See click documentation at http://click.pocoo.org/")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/climate_toolbox/climate_toolbox.py b/climate_toolbox/climate_toolbox.py
index 322099a..d843cbb 100644
--- a/climate_toolbox/climate_toolbox.py
+++ b/climate_toolbox/climate_toolbox.py
@@ -305,3 +305,240 @@ def _prepare_spatial_weights_data(weights_file=None):
         inplace=True)
 
     return df
+
+
+def _reindex_spatial_data_to_regions(ds, df):
+    '''
+    Reindexes spatial and segment weight data to regions
+    Enables region index-based math operations
+    Parameters
+    ----------
+    ds: xarray Dataset
+    df: pandas DataFrame
+    Returns
+    -------
+    Xarray DataArray
+    '''
+
+    # use vectorized indexing in xarray >= 0.10
+    if LooseVersion(xr.__version__) > LooseVersion('0.9.999'):
+
+        lon_indexer = xr.DataArray(df.lon.values, dims=('reshape_index', ))
+        lat_indexer = xr.DataArray(df.lat.values, dims=('reshape_index', ))
+
+        return ds.sel(lon=lon_indexer, lat=lat_indexer)
+
+    else:
+        res = ds.sel_points(
+            'reshape_index',
+            lat=df.lat.values,
+            lon=df.lon.values)
+
+        return res
+
+
+def _aggregate_reindexed_data_to_regions(
+        ds,
+        variable,
+        aggwt,
+        agglev,
+        weights,
+        backup_aggwt='areawt'):
+    '''
+    Performs weighted avg for climate variable by region
+
+    Parameters
+    ----------
+
+    ds: xarray.DataArray
+
+    variable: str
+        name of the data variable
+
+    aggwt: str
+        variable to weight by (i.e popwt, areawt, cropwt)
+
+    agglev: str
+        indicates which regional id scheme to select in the dataframe
+
+    weight: pd.DataFrame
+        pandas DataFrame of weights
+
+    backup_aggwt: str, optional
+        aggregation weight to use in regions with no aggwt data (default
+        'areawt')
+
+    '''
+
+    ds.coords[agglev] = xr.DataArray(
+                weights[agglev].values,
+                dims={'reshape_index': weights.index.values})
+
+    # format weights
+    ds[aggwt] = xr.DataArray(
+                weights[aggwt].values,
+                dims={'reshape_index': weights.index.values})
+
+    ds[aggwt] = (
+        ds[aggwt]
+        .where(ds[aggwt] > 0)
+        .fillna(weights[backup_aggwt].values))
+
+    weighted = xr.Dataset({
+        variable: (
+            (
+                (ds[variable]*ds[aggwt])
+                .groupby(agglev)
+                .sum(dim='reshape_index')) /
+            (
+                ds[aggwt]
+                .groupby(agglev)
+                .sum(dim='reshape_index')))})
+
+    return weighted
+
+
+'''
+================
+Public Functions
+================
+'''
+
+
+def load_bcsd(fp, varname, lon_name='lon', broadcast_dims=('time',)):
+    '''
+    Read and prepare climate data
+
+    After reading data, this method also fills NA values using linear
+    interpolation, and standardizes longitude to -180:180
+
+    Parameters
+    ----------
+    fp: str
+        File path or dataset
+
+    varname: str
+        Variable name to be read
+
+    lon_name : str, optional
+        Name of the longitude dimension (defualt selects from ['lon' or
+        'longitude'])
+
+    Returns
+    -------
+    xr.Dataset
+         xarray dataset loaded into memory
+    '''
+
+    if lon_name is not None:
+        lon_names = [lon_name]
+
+    if hasattr(fp, 'sel_points'):
+        ds = fp
+
+    else:
+        with xr.open_dataset(fp) as ds:
+            ds.load()
+
+    _fill_holes_xr(ds, varname, broadcast_dims=broadcast_dims)
+    return _standardize_longitude_dimension(ds, lon_names=lon_names)
+
+
+def load_baseline(fp, varname, lon_name='lon', broadcast_dims=None):
+    '''
+    Read and prepare climate data
+
+    After reading data, this method also fills NA values using linear
+    interpolation, and standardizes longitude to -180:180
+
+    Parameters
+    ----------
+    fp: str
+        File path or dataset
+
+    varname: str
+        Variable name to be read
+
+    lon_name : str, optional
+        Name of the longitude dimension (defualt selects from ['lon' or
+        'longitude'])
+
+    Returns
+    -------
+    xr.Dataset
+         xarray dataset loaded into memory
+    '''
+
+    if lon_name is not None:
+        lon_names = [lon_name]
+
+    if broadcast_dims is None:
+        broadcast_dims = tuple([])
+
+    if hasattr(fp, 'sel_points'):
+        ds = fp
+
+    else:
+        with xr.open_dataset(fp) as ds:
+            ds.load()
+
+    if 'lat' in ds.data_vars:
+        ds = ds.set_coords('lat')
+        ds = ds.swap_dims({'nlat': 'lat'})
+
+    if 'lon' in ds.data_vars:
+        ds = ds.set_coords('lon')
+        ds = ds.swap_dims({'nlon': 'lon'})
+
+    _fill_holes_xr(ds, varname, broadcast_dims=broadcast_dims)
+    return _standardize_longitude_dimension(ds, lon_names=lon_names)
+
+
+def weighted_aggregate_grid_to_regions(
+        ds,
+        variable,
+        aggwt,
+        agglev,
+        weights=None):
+    '''
+    Computes the weighted reshape of gridded data
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        xarray Dataset to be aggregated. Must have 'lat' and 'lon' in the
+        coordinates.
+
+    variable : str
+        name of the variable to be aggregated
+
+    aggwt : str
+        Weighting variable (e.g. 'popwt', 'areawt'). This must be a column name
+        in the weights file.
+
+    agglev : str
+        Target regional aggregation level (e.g. 'ISO', 'hierid'). This must be
+        a column name in the weights file.
+
+    weights : str, optional
+        Regional aggregation weights (default agglomerated-world-new BCSD
+        segment weights)
+
+    Returns
+    -------
+    ds: xr.Dataset
+        weighted and averaged dataset based on agglev
+    '''
+
+    if weights is None:
+        weights = _prepare_spatial_weights_data()
+
+    ds = _reindex_spatial_data_to_regions(ds, weights)
+    ds = _aggregate_reindexed_data_to_regions(
+        ds,
+        variable,
+        aggwt,
+        agglev,
+        weights)
+
+    return ds
diff --git a/climate_toolbox/aggregations/__init__.py b/climate_toolbox/geo/__init__.py
similarity index 100%
rename from climate_toolbox/aggregations/__init__.py
rename to climate_toolbox/geo/__init__.py
diff --git a/climate_toolbox/geo/distance.py b/climate_toolbox/geo/distance.py
new file mode 100644
index 0000000..494fd96
--- /dev/null
+++ b/climate_toolbox/geo/distance.py
@@ -0,0 +1,64 @@
+
+#             model             major (km)   minor (km)     flattening
+ELLIPSOIDS = {'WGS-84':        (6378.137,    6356.7523142,  1 / 298.257223563),
+              'GRS-80':        (6378.137,    6356.7523141,  1 / 298.257222101),
+              'Airy (1830)':   (6377.563396, 6356.256909,   1 / 299.3249646),
+              'Intl 1924':     (6378.388,    6356.911946,   1 / 297.0),
+              'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465),
+              'GRS-67':        (6378.1600,   6356.774719,   1 / 298.25),
+              }
+
+
+EARTH_RADIUS = 6371.009
+
+def great_circle(ax, ay, bx, by, radius=EARTH_RADIUS):
+    '''
+    calculate the great circle distance (km) between points
+
+    Provide points (ax, ay) and (bx, by) as floats, or as
+    vectors. If ax and ay are vectors or arrays of the
+    same shape, the element-wise distance will be found
+    between points in the vectors/arrays. If ax, ay are
+    (Mx1) column vectors and (bx, by) are (1xN) row
+    vectors, the vectors will be broadcast using numpy
+    broadcasting rules and the distance between each pair
+    of points will be returned as an (MxN) matrix.
+
+    Parameters
+    -----------
+    ax : float or array
+        x/long of point a
+    ay : float or array
+        y/lat of point a
+    bx : float or array
+        x/long of point b
+    by : float or array
+        y/lat of point b
+    radius : float, optional
+        Radius of the sphere on which to calculate the great
+        circle distance (default is to use the Earth's radius in
+        km, `6371.009`). Values returned will be in units of the
+        radius provided.
+
+    Returns
+    --------
+    distance : float or array
+        great circle distance between points a and b. Units will
+        match the radius provided (default km)
+    '''
+
+    lat1, lng1 = np.radians(ay), np.radians(ax)
+    lat2, lng2 = np.radians(by), np.radians(bx)
+
+    sin_lat1, cos_lat1 = np.sin(lat1), np.cos(lat1)
+    sin_lat2, cos_lat2 = np.sin(lat2), np.cos(lat2)
+
+    delta_lng = lng2 - lng1
+    cos_delta_lng, sin_delta_lng = np.cos(delta_lng), np.sin(delta_lng)
+
+    d = np.arctan2(np.sqrt((cos_lat2 * sin_delta_lng) ** 2 +
+                   (cos_lat1 * sin_lat2 -
+                    sin_lat1 * cos_lat2 * cos_delta_lng) ** 2),
+              sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)
+
+    return radius * d
diff --git a/climate_toolbox/io/__init__.py b/climate_toolbox/io/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/climate_toolbox/io/io.py b/climate_toolbox/io/io.py
deleted file mode 100644
index 0043d74..0000000
--- a/climate_toolbox/io/io.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import xarray as xr
-
-def load_bcsd(fp, varname, lon_name='lon', broadcast_dims=('time',)):
-    '''
-    Read and prepare climate data
-
-    After reading data, this method also fills NA values using linear
-    interpolation, and standardizes longitude to -180:180
-
-    Parameters
-    ----------
-    fp: str
-        File path or dataset
-
-    varname: str
-        Variable name to be read
-
-    lon_name : str, optional
-        Name of the longitude dimension (defualt selects from ['lon' or
-        'longitude'])
-
-    Returns
-    -------
-    xr.Dataset
-         xarray dataset loaded into memory
-    '''
-
-    if lon_name is not None:
-        lon_names = [lon_name]
-
-    if hasattr(fp, 'sel_points'):
-        ds = fp
-
-    else:
-        with xr.open_dataset(fp) as ds:
-            ds.load()
-
-    _fill_holes_xr(ds, varname, broadcast_dims=broadcast_dims)
-    return _standardize_longitude_dimension(ds, lon_names=lon_names)
-
-
-def load_gmfd(fp, varname, lon_name='lon', broadcast_dims=('time',)):
-    pass
-
-
-def load_best(fp, varname, lon_name='lon', broadcast_dims=('time',)):
-    pass
diff --git a/climate_toolbox/transformations/__init__.py b/climate_toolbox/transformations/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/climate_toolbox/transformations/transformations.py b/climate_toolbox/transformations/transformations.py
deleted file mode 100644
index 18de472..0000000
--- a/climate_toolbox/transformations/transformations.py
+++ /dev/null
@@ -1,102 +0,0 @@
-import xarray as xr
-import numpy as np
-
-from climate_toolbox.utils.utils import remove_leap_days
-
-def edd_ag(ds_tasmax, ds_tasmin, threshold):
-    '''
-    Note: there are implicitly three cases:
-
-        1. tmax > threshold & tmin < threshold
-        2. tmax > threshold & tmin > threshold
-        3. tmax <= threshold
-    
-    Case (1) is the first part of the np.where() statement.
-    Case (2) is also the first part of this statement, which returns 'NA' in
-    this case and so is not included in the final summation of HDD. (I.e., in
-    case (2), the HDD should be zero. Instead we get 'NA', which functions as
-    the equivalent of a zero value in the subsequent code.)
-    Case (3) is, of course, the second part of the np.where() statement.
-    
-    Parameters
-    ----------
-    ds : Dataset
-        xarray.Dataset with two variables: tasmin and tasmax. 
-        tasmin and tasmax are in Kelvin and are indexed by
-        impact region (``hierid``) and day (``time``) in a
-        365-day calendar.
-        
-    Returns
-    -------
-    ds : Dataset
-        xarray.Dataset with dimensions ``(hierid, threshold)``
-    '''
-
-    # convert from K to C
-    tmax = (ds_tasmax.tasmax - 273.15)
-    tmin = (ds_tasmin.tasmin - 273.15)
-
-    #get the 
-    snyder_m = (tmax + tmin)/2
-    snyder_w = (tmax - tmin)/2
-    snyder_theta = np.arcsin( (threshold - snyder_m)/snyder_w )
-
-    transdata = np.where(
-        tmin.values < threshold,
-            np.where(
-                tmax.values > threshold,
-            ((snyder_m.values - threshold) * (np.pi/2 - snyder_theta.values) +
-                         snyder_w.values * np.cos(snyder_theta.values) ) / np.pi, 0),
-        snyder_m.values - threshold)
-
-    
-    return xr.DataArray(transdata, dims=tmax.dims, coords=tmax.coords)
-
-
-def tas_poly(ds, power):
-    '''
-    Daily average temperature (degrees C), raised to a power
-
-    Leap years are removed before counting days (uses a 365 day
-    calendar).
-    '''
-
-    powername = ordinal(power)
-
-    description = format_docstr(('''
-            Daily average temperature (degrees C){raised}
-
-            Leap years are removed before counting days (uses a 365 day
-            calendar).
-            '''.format(
-                raised='' if power == 1 else (
-                    ' raised to the {powername} power'
-                    .format(powername=powername)))).strip())
-
-    ds1 = xr.Dataset()
-
-    # remove leap years
-    ds = remove_leap_days(ds)
-
-    # do transformation
-    ds1[varname] = (ds.tas - 273.15)**power
-
-    # Replace datetime64[ns] 'time' with YYYYDDD int 'day'
-    if ds.dims['time'] > 365:
-        raise ValueError
-
-    ds1.coords['day'] = ds['time.year']*1000 + np.arange(1, len(ds.time)+1)
-    ds1 = ds1.swap_dims({'time': 'day'})
-    ds1 = ds1.drop('time')
-
-    ds1 = ds1.rename({'day': 'time'})
-
-    # document variable
-    ds1[varname].attrs['units'] = (
-        'C^{}'.format(power) if power > 1 else 'C')
-
-    ds1[varname].attrs['long_title'] = description.splitlines()[0]
-    ds1[varname].attrs['description'] = description
-    ds1[varname].attrs['variable'] = varname
-
-    return ds1
\ No newline at end of file
diff --git a/climate_toolbox/utils/__init__.py b/climate_toolbox/utils/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/climate_toolbox/utils/utils.py b/climate_toolbox/utils/utils.py
deleted file mode 100644
index 59d1b97..0000000
--- a/climate_toolbox/utils/utils.py
+++ /dev/null
@@ -1,195 +0,0 @@
-'''
-Handy functions for standardizing the format of climate data
-'''
-
-import xarray as xr
-import numpy as np
-
-
-def convert_lons_mono(ds, lon_name='longitude'):
-    ''' Convert longitude from -180-180 to 0-360 '''
-    ds[lon_name].values = np.where(
-        ds[lon_name].values < 0, 360 + ds[lon_name].values, ds[lon_name].values
-        )
-
-    # sort the dataset by the new lon values
-    ds = ds.sel(**{lon_name: np.sort(ds[lon_name].values)})
-
-    return ds
-
-
-def convert_lons_split(ds, lon_name='longitude'):
-    ''' Convert longitude from 0-360 to -180-180 '''
-    ds[lon_name].values = xr.where(
-        ds[lon_name] > 180, ds[lon_name] - 360, ds[lon_name])
-
-    # sort the dataset by the new lon values
-    ds = ds.sel(**{lon_name: np.sort(ds[lon_name].values)})
-
-    return ds
-
-
-def rename_coords_to_lon_and_lat(ds):
-    ''' Rename Dataset spatial coord names to:
-        lat, lon
-    '''
-    if 'latitude' in ds.coords:
-        ds = ds.rename({'latitude': 'lat'})
-    if 'longitude' in ds.coords:
-        ds = ds.rename({'longitude': 'lon'})
-    elif 'long' in ds.coords:
-        ds = ds.rename({'long': 'lon'})
-
-    if 'z' in ds.coords:
-        ds = ds.drop('z').squeeze()
-
-    return ds
-
-
-def rename_coords_to_longitude_and_latitude(ds):
-    ''' Rename Dataset spatial coord names to:
-        latitude, longitude
-    '''
-    if 'lat' in ds.coords:
-        ds = ds.rename({'lat': 'latitude'})
-    if 'lon' in ds.coords:
-        ds = ds.rename({'lon': 'longitude'})
-    elif 'long' in ds.coords:
-        ds = ds.rename({'long': 'longitude'})
-
-    if 'z' in ds.coords:
-        ds = ds.drop('z').squeeze()
-
-    return ds
-
-
-def remove_leap_days(ds):
-    ds = ds.loc[{
-        'time': ~((ds['time.month'] == 2) & (ds['time.day'] == 29))}]
-
-    return ds
-
-
-def season_boundaries(growing_days):
-    ''' Returns the sorted start and end date of growing season
-    '''
-
-    # the longitude values of the data is off, we need to scale it
-    growing_days.longitude.values = growing_days.longitude.values - 180
-    # we then sort by longitude
-    growing_days = growing_days.sortby('longitude')
-
-    # construct the ds
-    gdd_sorted = xr.DataArray(
-        # xarray has no method to sort along an axis
-        # we use np.sort but construct the matrix from a xarray dataArray
-        # we use transpose to track the axis we want to sort along
-        np.sort(
-            growing_days.variable.transpose(
-                'latitude', 'longitude', 'z').values, axis=2),
-        dims=('latitude', 'longitude', 'sort'),
-        coords={
-                'latitude': growing_days.latitude,
-                'longitude': growing_days.longitude,
-                'sort': pd.Index(['min', 'max'])
-                }
-        )
-
-    # we can then select an axis in the sorted dataarray as min
-    min_day, max_day = gdd_sorted.sel(sort='min'), gdd_sorted.sel(sort='max')
-
-    return min_day, max_day
-
-
-def get_daily_growing_season_mask(lat, lon, time, growing_days_path):
-    '''
-
-    Constructs a mask for days in the within calendar growing season
-
-    Parameters
-    ----------
-    lat: xr.DataArray coords object
-    lon: xr.DataArray coords object
-    time: xr.DataArray coords object
-    growing_days_path: str
-
-    Returns
-    -------
-    DataArray
-        xr.DataArray of masked lat x lon x time
-
-    '''
-
-    growing_days = xr.open_dataset(growing_days_path)
-
-    # find the min and max for the growing season
-    min_day, max_day = season_boundaries(growing_days)
-
-    data = np.ones((lat.shape[0], lon.shape[0], time.shape[0]))
-    # create an array of ones in the shape of the data
-    ones = xr.DataArray(
-                    data, coords=[lat, lon, time], dims=['lat', 'lon', 'time'])
-
-    # mask the array around the within calendar year start and end times
-    # of growing season
-    mask = (
-        (ones['time.dayofyear'] >= min_day) &
-        (ones['time.dayofyear'] <= max_day))
-
-    # apply this mask and
-    finalmask = (
-        mask.where(
-            growing_days.variable.sel(z=2) >=
-            growing_days.variable.sel(z=1)).fillna(1-mask).where(
-                ~growing_days.variable.sel(z=1, drop=True).isnull()
-                ).rename({'latitude': 'lat', 'longitude': 'lon'})
-            )
-
-    return finalmask
-
-
-def edd_ag(ds_tasmax, ds_tasmin, threshold):
-    '''
-
-    Note: there are implicitly three cases:
-
-        1. tmax > threshold & tmin < threshold
-        2. tmax > threshold & tmin > threshold
-        3. tmax <= threshold
-
-    Case (1) is the first part of the np.where() statement.
-    Case (2) is also the first part of this statement, which returns 'NA' in
-    this case and so is not included in the final summation of HDD. (I.e., in
-    case (2), the HDD should be zero. Instead we get 'NA', which functions as
-    the equivalent of a zero value in the subsequent code.)
-    Case (3) is, of course, the second part of the np.where() statement.
-
-    Parameters
-    ----------
-    ds : Dataset
-        xarray.Dataset with two variables: tasmin and tasmax.
-        tasmin and tasmax are in Kelvin and are indexed by
-        impact region (``hierid``) and day (``time``) in a
-        365-day calendar.
-
-    Returns
-    -------
-    ds : Dataset
-        xarray.Dataset with dimensions ``(hierid, threshold)``
-    '''
-
-    # convert from K to C
-    tmax = (ds_tasmax.tasmax - 273.15)
-    tmin = (ds_tasmin.tasmin - 273.15)
-
-    snyder_m = (tmax + tmin)/2
-    snyder_w = (tmax - tmin)/2
-    snyder_theta = np.arcsin((threshold - snyder_m)/snyder_w)
-
-    transdata = np.where(
-        tmin.values < threshold, np.where(tmax.values > threshold, (
-            (snyder_m.values - threshold) * (np.pi/2 - snyder_theta.values) +
-            snyder_w.values * np.cos(snyder_theta.values)
-            ) / np.pi, 0), snyder_m.values - threshold)
-
-    return xr.DataArray(transdata, dims=tmax.dims, coords=tmax.coords)
diff --git a/docs/.gitignore b/docs/.gitignore
deleted file mode 100644
index 66385f4..0000000
--- a/docs/.gitignore
+++ /dev/null
@@ -1,3 +0,0 @@
-/climate_toolbox.rst
-/climate_toolbox.*.rst
-/modules.rst
diff --git a/docs/climate_toolbox.geo.rst b/docs/climate_toolbox.geo.rst
new file mode 100644
index 0000000..1b46656
--- /dev/null
+++ b/docs/climate_toolbox.geo.rst
@@ -0,0 +1,22 @@
+climate\_toolbox.geo package
+============================
+
+Submodules
+----------
+
+climate\_toolbox.geo.distance module
+------------------------------------
+
+.. automodule:: climate_toolbox.geo.distance
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: climate_toolbox.geo
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/docs/climate_toolbox.rst b/docs/climate_toolbox.rst
new file mode 100644
index 0000000..e966455
--- /dev/null
+++ b/docs/climate_toolbox.rst
@@ -0,0 +1,37 @@
+climate\_toolbox package
+========================
+
+Subpackages
+-----------
+
+.. toctree::
+
+    climate_toolbox.geo
+
+Submodules
+----------
+
+climate\_toolbox.cli module
+---------------------------
+
+.. automodule:: climate_toolbox.cli
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+climate\_toolbox.climate\_toolbox module
+----------------------------------------
+
+.. automodule:: climate_toolbox.climate_toolbox
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Module contents
+---------------
+
+.. automodule:: climate_toolbox
+    :members:
+    :undoc-members:
+    :show-inheritance:
diff --git a/docs/modules.rst b/docs/modules.rst
new file mode 100644
index 0000000..ece927a
--- /dev/null
+++ b/docs/modules.rst
@@ -0,0 +1,7 @@
+climate_toolbox
+===============
+
+.. toctree::
+   :maxdepth: 4
+
+   climate_toolbox
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index eec9209..0000000
--- a/requirements.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-pytest-cov==2.7.1
-scipy==1.3.0
-xarray==0.12.1
-toolz==0.9.0
diff --git a/requirements_dev.txt b/requirements_dev.txt
index b3b3a96..809081b 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -11,3 +11,4 @@ pytest==3.1.3
 pytest-runner==2.11.1
 pytest-cov==2.5.1
 py==1.4.34
+


Test patch:
diff --git a/tests/test_climate_toolbox.py b/tests/test_climate_toolbox.py
index fa2dd35..fe8b614 100644
--- a/tests/test_climate_toolbox.py
+++ b/tests/test_climate_toolbox.py
@@ -6,8 +6,6 @@
 import pytest
 
 from climate_toolbox import climate_toolbox as ctb
-from climate_toolbox.utils.utils import *
-from climate_toolbox.aggregations.aggregations import _reindex_spatial_data_to_regions, _aggregate_reindexed_data_to_regions
 
 import numpy as np
 import pandas as pd
@@ -108,7 +106,7 @@ def test_reindex_spatial_weights(clim_data, weights):
 
     assert not clim_data.temperature.isnull().any()
 
-    ds = _reindex_spatial_data_to_regions(clim_data, weights)
+    ds = ctb._reindex_spatial_data_to_regions(clim_data, weights)
 
     assert ds.temperature.shape == (len(ds['lon']), len(ds['time']))
     assert 'reshape_index' in ds.dims
@@ -117,94 +115,15 @@ def test_reindex_spatial_weights(clim_data, weights):
 def test_weighting(clim_data, weights):
 
     assert np.isnan(weights['popwt'].values).any()
-    ds = _reindex_spatial_data_to_regions(clim_data, weights)
+    ds = ctb._reindex_spatial_data_to_regions(clim_data, weights)
     assert not ds.temperature.isnull().any()
 
-    wtd = _aggregate_reindexed_data_to_regions(
+    wtd = ctb._aggregate_reindexed_data_to_regions(
         ds, 'temperature', 'popwt', 'ISO', weights)
 
     assert not wtd.temperature.isnull().any()
 
-    wtd = _aggregate_reindexed_data_to_regions(
+    wtd = ctb._aggregate_reindexed_data_to_regions(
         ds, 'temperature', 'areawt', 'ISO', weights)
 
     assert not wtd.temperature.isnull().any()
-
-
-def test_rename_coords_to_lon_and_lat():
-    ds = xr.Dataset(
-        coords={'z': [1.20, 2.58], 'long': [156.6, 38.48]})
-
-    ds = rename_coords_to_lon_and_lat(ds)
-    coords = ds.coords
-
-    assert 'z' not in coords
-    assert coords.z is None
-    assert 'lon' in coords and 'long' not in coords
-
-
-def test_rename_coords_to_lon_and_lat():
-    ds = xr.Dataset(
-        coords={'latitude': [71.32, 72.58], 'longitude': [156.6, 38.48]})
-
-    ds = rename_coords_to_lon_and_lat(ds)
-    coords = ds.coords
-
-    assert 'lat' in coords and 'latitude' not in coords
-    assert 'lon' in coords and 'longitude' not in coords
-
-
-def test_rename_coords_to_longitude_and_latitude():
-    ds = xr.Dataset(
-        coords={'lat': [71.32, 72.58], 'lon': [156.6, 38.48]})
-    ds = rename_coords_to_longitude_and_latitude(ds)
-    coords = ds.coords
-
-    assert 'latitude' in coords and 'lat' not in coords
-    assert 'longitude' in coords and 'lon' not in coords
-
-
-def test_rename_coords_to_longitude_and_latitude_with_clim_data(clim_data):
-    ds = rename_coords_to_longitude_and_latitude(clim_data)
-    coords = ds.coords
-
-    assert 'latitude' in coords and 'lat' not in coords
-    assert 'longitude' in coords and 'lon' not in coords
-
-
-def test_convert_lons_mono():
-    ds = xr.Dataset(coords={'lon': [-156.6, -38.48]})
-    expected = np.array([203.4, 321.52])
-
-    ds = convert_lons_mono(ds, lon_name='lon')
-
-    np.testing.assert_array_equal(ds.lon.values, expected)
-
-
-def test_convert_lons_split():
-    ds = xr.Dataset(coords={'longitude': [300, 320]})
-    expected = np.array([-60, -40])
-
-    ds = convert_lons_split(ds)
-
-    np.testing.assert_array_equal(ds.longitude.values, expected)
-
-
-def test_remove_leap_days():
-    da = xr.DataArray(
-        np.random.rand(4, 3),
-        [('time', pd.date_range('2000-02-27', periods=4)),
-         ('space', ['IA', 'IL', 'IN'])])
-    leap_day = np.datetime64('2000-02-29')
-
-    da = remove_leap_days(da)
-
-    assert leap_day not in da.coords['time'].values
-
-
-def test_remove_leap_days_with_clim_data(clim_data):
-    leap_day = np.datetime64('2000-02-29')
-
-    da = remove_leap_days(clim_data)
-
-    assert leap_day not in da.coords['time'].values