Repository: CSIRO-Precision-Agriculture/pyprecag (language: python)

Issue:
pyprecag for Python 3
The migration of pyprecag to support python 3 is underway and expected to be released in July 2020
Geopandas 0.6.0 causing TypeErrors
There are TypeErrors in describe.py, causing tests to error and fail. 
The error was as follows:
```
File "pyprecag/describe.py", line 318, in get_column_properties
fldtype = type(np.asscalar(np.zeros(1, _type))).name
TypeError: data type not understood
``'


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 740d9ef..56140cd 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,5 +1,10 @@
 Changelog
 =========
+** 1.0.0 (2020-08-07)**
+ * #40 Port to Python3
+ * Updated Blockgrid to support a group-by batch mode.
+ * ``block_size`` argument removed from ``kriging_ops.prepare_for_vesper_krige``
+
 **0.3.1 (2020-05-27)**
  Bug Fixes
   * Support for Geopandas (upto 0.5.1)
diff --git a/CITATION.rst b/CITATION.rst
index 768a2e3..5c6fbb9 100644
--- a/CITATION.rst
+++ b/CITATION.rst
@@ -1,4 +1,4 @@
 If you use pyprecag for any published work, please cite it using the reference below:
 
 Ratcliff, Christina; Gobbett, David; Bramley, Rob (2019): pyprecag  - a Python package for the analysis of
-precision agriculture data.. v1. CSIRO. Software Collection. https://doi.org/10.25919/5c731a41954ce
+precision agriculture data.. v1. CSIRO. Software Collection. https://doi.org/10.25919/5d23f6ba8cff0
diff --git a/LICENSE b/LICENSE
index d627e9d..eed7e3e 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 CSIRO Open Source Software License Agreement (variation of the BSD / MIT License)
 
-Copyright (c) 2018, Commonwealth Scientific and Industrial Research Organisation (CSIRO) ABN 41 687 119 230.
+Copyright (c) 2018-2020, Commonwealth Scientific and Industrial Research Organisation (CSIRO) ABN 41 687 119 230.
 
 All rights reserved. CSIRO is willing to grant you a license to this pyprecag on the following terms, except where otherwise indicated for third party material.
 
diff --git a/README.rst b/README.rst
index ea57a12..cadd4e1 100644
--- a/README.rst
+++ b/README.rst
@@ -13,17 +13,18 @@ A suite of tools for Precision Agriculture data analysis
 * Version: 0.3.1
 * License: `CSIRO Open Source Software License Agreement <https://github.com/CSIRO-Precision-Agriculture/pyprecag/blob/master/LICENSE>`__
 
-pyprecag is currently supported only on Python version 2.7
+pyprecag is currently supported for Python version 2.7 and 3.7
 
 .. include:: ./CITATION.rst
 
 Dependencies
 ------------
 
-* Python version 2.7
+* Python version 2.7, 3.7
 * `GDAL <https://www.gdal.org/>`_
 * `Fiona <https://github.com/Toblerity/Fiona>`_
 * `Rasterio <https://github.com/mapbox/rasterio>`_
+* `Geopandas <https://github.com/geopandas/geopandas>`_
 * `VESPER <https://sydney.edu.au/agriculture/pal/software/vesper.shtml>`_ (for Kriging)
 
 .. include:: ./docs/installation.rst
diff --git a/docs/conf.py b/docs/conf.py
index 4fa424f..99c7b25 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -28,7 +28,7 @@
 source_suffix = '.rst'
 master_doc = 'index'
 project = 'pyprecag'
-copyright = '2018, CSIRO'
+copyright = '2018-2020, CSIRO'
 version = release = '0.3.1'
 pygments_style = 'perldoc'
 templates_path = ['.']
diff --git a/pyprecag/__init__.py b/pyprecag/__init__.py
index 7f4ddf6..b8a95b2 100644
--- a/pyprecag/__init__.py
+++ b/pyprecag/__init__.py
@@ -9,12 +9,17 @@
 import sys
 import subprocess
 import tempfile
+
+import six
+
 from . import config
 
 __author__ = 'Christina Ratcliff',
 __email__ = 'Christina.Ratcliff@csiro.au',
 __version__ = '0.3.1'
 
+number_types = six.integer_types + (float, )
+
 TEMPDIR = os.path.join(tempfile.gettempdir(), 'PrecisionAg')
 
 if not os.path.exists(TEMPDIR):
@@ -31,6 +36,11 @@
         process = subprocess.Popen(['gdal-config','--datadir'],
             stdout=subprocess.PIPE)
         output, error = process.communicate()
+        output = six.ensure_str(output)
+        # if there is no error in the subprocess, error is None
+        # which six refuses to coerce to str
+        if error is not None:
+            error = six.ensure_str(error)
         gdal_data = output.split('\n')[0]
         assert os.path.isfile(os.path.join(gdal_data, 'gcs.csv'))
     except:
diff --git a/pyprecag/bandops.py b/pyprecag/bandops.py
index e304c80..739afda 100644
--- a/pyprecag/bandops.py
+++ b/pyprecag/bandops.py
@@ -3,6 +3,7 @@
 import numpy as np
 from numpy import ma
 import rasterio
+import six
 
 
 class BandMapping(collections.MutableMapping, dict):
@@ -30,7 +31,9 @@ def __setitem__(self, key, value):
 
         if key not in allowed_keys:
             raise AttributeError(
-                'BandMapping has no attribute {}. Allowed keys are {}'.format(key, ', '.join(self.__defaults.keys())))
+                'BandMapping has no attribute {}. Allowed keys are {}'.format(
+                    key, ', '.join(sorted(self.__defaults.keys())))
+                )
 
         if not isinstance(value, int):
             raise ValueError('{v} is not an integer'.format(v=value))
@@ -118,11 +121,11 @@ def calculate(self, index_name, raster, src_nodata=None, dest_nodata=-9999):
             pass
         elif not os.path.exists(raster):
             raise ValueError('Image file does not exist')
-
-        if src_nodata is not None and not isinstance(src_nodata, (int, float, long)):
+        number_types = six.integer_types + (float,)
+        if src_nodata is not None and not isinstance(src_nodata, number_types):
             raise ValueError('src_nodata should be numeric (int, float or long)')
 
-        if not dest_nodata and not isinstance(dest_nodata, (int, float, long)):
+        if not dest_nodata and not isinstance(dest_nodata, number_types):
             raise ValueError('dst_nodata should be numeric (int, float or long)')
 
         self.__image = raster
diff --git a/pyprecag/config.json b/pyprecag/config.json
index 7f23cb9..f699686 100644
--- a/pyprecag/config.json
+++ b/pyprecag/config.json
@@ -1,6 +1,7 @@
 {
     "crsLookupURL": "http://prj2epsg.org/search.json",
     "geoCSV": {
+        "file_encoding": "utf-8",
         "xCoordinate_ColumnName": [
             "LONGITUDE",
             "LON",
diff --git a/pyprecag/convert.py b/pyprecag/convert.py
index c4141ad..0f9b593 100644
--- a/pyprecag/convert.py
+++ b/pyprecag/convert.py
@@ -14,6 +14,8 @@
 import pandas as pd
 from pandas.core.dtypes.common import is_string_dtype
 
+import six
+
 from osgeo import ogr, gdal
 import geopandas
 from geopandas import GeoDataFrame, GeoSeries
@@ -22,13 +24,13 @@
 from rasterio import features
 
 from fiona import collection as fionacoll
-from fiona.crs import from_epsg
 from shapely.geometry import Point, mapping, shape, LineString
 
 from . import crs as pyprecag_crs
+from pyprecag.crs import from_epsg
+
 from . import TEMPDIR, config
-from .describe import CsvDescribe, predictCoordinateColumnNames, VectorDescribe, \
-    save_geopandas_tofile
+from .describe import CsvDescribe, predictCoordinateColumnNames, VectorDescribe, save_geopandas_tofile
 
 from .errors import GeometryError
 from .raster_ops import raster_snap_extent, create_raster_transform
@@ -190,9 +192,9 @@ def convert_polygon_to_grid(in_shapefilename,
 
     # for pixel data types see http://www.gdal.org/gdal_8h.html#a22e22ce0a55036a96f652765793fb7a4
     tif_drv = gdal.GetDriverByName('GTiff')
-    if tif_drv is None:
-        # If the above Fails then here is an alternate version.
-        tif_drv = Find_GDALDriverByName('GTiff')
+    # if tif_drv is None:
+    #     # If the above Fails then here is an alternate version.
+    #     tif_drv = Find_GDALDriverByName('GTiff')
 
     raster_ds = tif_drv.Create(out_rasterfilename, x_cols, y_rows, 1, gdal.GDT_Int16)
 
@@ -218,8 +220,10 @@ def convert_polygon_to_grid(in_shapefilename,
         LOGGER.exception("Error rasterizing layer: {}".format(err), True)
 
     else:
-        LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                                   dur=timedelta(seconds=time.time() - start_time)))
+        LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(
+            inspect.currentframe().f_code.co_name, '',
+            dur=str(timedelta(seconds=time.time() - start_time))
+        ))
 
     # close the data sources
     r_band = None
@@ -261,8 +265,10 @@ def convert_grid_to_vesper(in_rasterfilename, out_vesperfilename):
                         xy = oRasterDS.xy(iRowY, iColX)
                         vesFile.write(' {}   {}\n'.format(xy[0], xy[1]))
 
-    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                               dur=timedelta(seconds=time.time() - start_time)))
+    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(
+        inspect.currentframe().f_code.co_name, '',
+        dur=str(timedelta(seconds=time.time() - start_time))
+    ))
     return
 
 
@@ -290,7 +296,7 @@ def add_point_geometry_to_dataframe(in_dataframe, coord_columns=None,
         raise TypeError("Input points should be a geopandas dataframe")
 
     for argCheck in [('coord_columns_epsg', coord_columns_epsg), ('out_epsg', out_epsg)]:
-        if not isinstance(argCheck[1], (int, long)):
+        if not isinstance(argCheck[1], six.integer_types):
             raise TypeError('{} must be a integer.'.format(argCheck[0]))
 
     if coord_columns is None:
@@ -312,15 +318,17 @@ def add_point_geometry_to_dataframe(in_dataframe, coord_columns=None,
         # insert feature id as column as first column and populate it using row number.
         in_dataframe.insert(0, 'FID', in_dataframe.index)
 
-    # combine lat and lon column to a shapely Point() object
-    in_dataframe['geometry'] = in_dataframe.apply(
-        lambda x: Point((float(x[x_column]), float(x[y_column]))), axis=1)
-
-    # Now, convert the pandas DataFrame into a GeoDataFrame. The geopandas constructor expects a
-    # geometry column which can consist of shapely geometry objects, so the column we created is
-    # just fine:
-    gdf_csv = geopandas.GeoDataFrame(in_dataframe, geometry='geometry',
-                                     crs=from_epsg(coord_columns_epsg))
+    if hasattr(geopandas,'points_from_xy'):   # implemented in geopandas 0.5.0
+        gdf_csv = geopandas.GeoDataFrame(in_dataframe,
+                                         geometry=geopandas.points_from_xy(in_dataframe[x_column],
+                                                                           in_dataframe[y_column]),
+                                         crs=from_epsg(coord_columns_epsg))
+    else:
+        # Convert the pandas DataFrame into a GeoDataFrame.
+        gdf_csv = geopandas.GeoDataFrame(in_dataframe,
+                                         geometry=[Point(x, y) for x, y in zip(in_dataframe[x_column],
+                                                                               in_dataframe[y_column])],
+                                         crs=from_epsg(coord_columns_epsg))
 
     # drop the original geometry columns to avoid confusion
     gdf_csv.drop([x_column, y_column], axis=1, inplace=True)
@@ -368,7 +376,7 @@ def convert_csv_to_points(in_csvfilename, out_shapefilename=None, coord_columns=
         raise IOError("Invalid path: {}".format(in_csvfilename))
 
     for argCheck in [('coord_columns_epsg', coord_columns_epsg), ('out_epsg', out_epsg)]:
-        if not isinstance(argCheck[1], (int, long)):
+        if not isinstance(argCheck[1], six.integer_types):
             raise TypeError('{} must be a integer.'.format(argCheck[0]))
 
     desc_csv = CsvDescribe(in_csvfilename)
@@ -385,8 +393,10 @@ def convert_csv_to_points(in_csvfilename, out_shapefilename=None, coord_columns=
 
         save_geopandas_tofile(gdf_csv, out_shapefilename, overwrite=True)
 
-    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                               dur=timedelta(seconds=time.time() - start_time)))
+    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(
+        str(inspect.currentframe().f_code.co_name), '',
+        dur=str(timedelta(seconds=time.time() - start_time))
+    ))
 
     return gdf_csv, gdf_crs
 
@@ -415,16 +425,15 @@ def convert_polygon_feature_to_raster(feature, pixel_size, value=1, nodata_val=-
     if not isinstance(feature, (pd.Series, geopandas.GeoSeries)):
         raise TypeError("Input feature should be a geopandas series")
 
-    if not isinstance(pixel_size, (int, long, float)):
+    if not isinstance(pixel_size, six.integer_types + (float, )):
         raise TypeError('Pixel size must be an integer or floating number.')
 
     if isinstance(value, str) and 'value' not in list(feature.index):
         raise TypeError('Value string {} is not a column name')
-    elif not isinstance(value, (int, long, float)):
+    elif not isinstance(value, six.integer_types + (float, )):
         raise TypeError('Value should be a column name, or a number')
 
-    transform, width, height, bbox = create_raster_transform(feature.geometry.bounds,
-                                                             pixel_size, snap_extent_to_pixel)
+    transform, width, height, bbox = create_raster_transform(feature.geometry.bounds, pixel_size, snap_extent_to_pixel)
 
     gdf_feat = GeoDataFrame(GeoSeries(feature.geometry),
                             geometry=GeoSeries(feature.geometry).geometry)
@@ -438,8 +447,7 @@ def convert_polygon_feature_to_raster(feature, pixel_size, value=1, nodata_val=-
     else:
         shapes = ((geom, value) for geom in gdf_feat.geometry)
 
-    burned = features.rasterize(shapes=shapes, out_shape=(height, width), fill=nodata_val,
-                                transform=transform)
+    burned = features.rasterize(shapes=shapes, out_shape=(height, width), fill=nodata_val, transform=transform)
     meta = {}
     meta.update(height=height, width=width,
                 transform=transform,
diff --git a/pyprecag/crs.py b/pyprecag/crs.py
index 2cf9172..e0c617b 100644
--- a/pyprecag/crs.py
+++ b/pyprecag/crs.py
@@ -2,14 +2,18 @@
 import logging
 
 import os
+
+import six
 import socket
 import warnings
-from urllib import urlencode
-import urllib2
 
-from fiona.crs import from_string, from_epsg
+from pkg_resources import parse_version
+from six.moves.urllib_parse import urlencode
+from six.moves.urllib_request import urlopen
+
 from osgeo import osr, gdal
 from shapely import geometry
+import pyproj
 
 from . import config
 from .errors import SpatialReferenceError
@@ -28,6 +32,7 @@ def __init__(self):
         self.epsg_number = None
         self.epsg = None
         self.epsg_predicted = False  # Did the epsg_number get set via the online lookup.
+        self.pyproj_version=pyproj.__version__
 
     def set_epsg(self, code):
         """Given an integer code, set the epsg_number and the EPSG-like mapping.
@@ -37,20 +42,20 @@ def set_epsg(self, code):
             code (int): An integer representing the EPSG code
 
          """
-        if code is None:
-            warnings.warn('EPSG Code is None and cant be set.')
+        if code is None or code == 0:
+            warnings.warn('EPSG Code is None or 0 and cant be set.')
             return
 
-        if isinstance(code, str):
+        if isinstance(code, six.string_types):
             try:
                 code = int(code)
             except:
-                raise TypeError('EPSG code must be a positive integer')
+                raise TypeError('EPSG code ({}) must be a positive integer'.format(code))
 
         if not isinstance(code, int):
-            raise TypeError('EPSG Code must be a positive integer')
+            raise TypeError('EPSG Code ({}) must be a positive integer'.format(code))
         if code <= 0:
-            raise ValueError('EPSG Code must be a positive integer')
+            raise ValueError('EPSG Code ({}) must be a positive integer above 0'.format(code))
 
         self.epsg = from_epsg(code)
         self.epsg_number = code
@@ -72,22 +77,25 @@ def getFromWKT(self, crs_wkt, bUpdateByEPSG=False):
 
         source = osr.SpatialReference()
         source.ImportFromWkt(crs_wkt)
+
         if source is not None:
             self.crs_wkt = crs_wkt
             self.srs = source
             self.proj4 = source.ExportToProj4()
 
             source.AutoIdentifyEPSG()
-
             code = source.GetAuthorityCode(None)
 
+            if code is None:
+                if parse_version(pyproj.__version__) >= parse_version('2.5.0'):
+                    pyproj_crs = pyproj.CRS(source.GetName())
+                    code = pyproj_crs.to_epsg()
             if code is None:
                 code = self.getEPSGFromSRS(source,  bUpdateToCorrectDefn=bUpdateByEPSG)
-            else:
-                code = int(source.GetAuthorityCode(None))
 
-            if code is not None:
-                self.set_epsg(code)
+            if code:
+                self.set_epsg(int(code))
+
 
     def getFromEPSG(self, epsg):
         """Create OGR Spatial Reference Object for an epsg_number number
@@ -98,13 +106,13 @@ def getFromEPSG(self, epsg):
             osgeo.osr.SpatialReference: The Spatial Reference Object for the input EPSG
 
         """
-        if isinstance(epsg, (str, unicode)):
+        if isinstance(epsg, six.string_types):
             try:
                 epsg = int(epsg.replace('EPSG:', ''))
             except:
                 raise TypeError('EPSG must be a number. Got - {}'.format(epsg))
 
-        if not isinstance(epsg, (int, long)):
+        if not isinstance(epsg, six.integer_types):
             raise TypeError('EPSG must be a number. Got - {}'.format(epsg))
 
         if epsg is None or epsg == 0:
@@ -143,7 +151,7 @@ def getEPSGFromSRS(self, osr_srs, bOnlineLookup=False, bUpdateToCorrectDefn=Fals
             int: EPSG Number for the matched Spatial Reference System.
 
         """
-        warnings.warn('bOnlineLookup is deprecated', PendingDeprecationWarning)
+        warnings.warn('bOnlineLookup is deprecated', DeprecationWarning)
 
         if osr_srs is None:
             return
@@ -183,7 +191,7 @@ def getEPSGFromSRS(self, osr_srs, bOnlineLookup=False, bUpdateToCorrectDefn=Fals
             import sqlite3
             import pandas as pd
 
-            connection =  sqlite3.connect(crs_database)
+            connection = sqlite3.connect(crs_database)
             cursor = connection.cursor()
             cs_df = pd.read_sql_query('Select  * from  alias_name', connection)
             del cursor
@@ -199,10 +207,31 @@ def getEPSGFromSRS(self, osr_srs, bOnlineLookup=False, bUpdateToCorrectDefn=Fals
 
         if bUpdateToCorrectDefn and epsg > 0:
             self.getFromEPSG(int(epsg))
+            self.set_epsg(epsg)
 
         return int(epsg)
 
 
+def from_epsg(epsg_number):
+    """
+    Get a coordinate system for the epsg number.
+    if pyproj is pre 2.5.0 then use a string created using fiona.crs.from_epsg, otherwise use the pyproj option.
+    Args:
+        epsg_number: The epsg number representing the coordinate system.
+
+    Returns:
+        object: pyproj.crs     The newer pyproj 2.5+ crs object compatible with Proj 6+
+                or
+               string          Derived from epsg using fiona.
+    """
+    if parse_version(pyproj.__version__) >= parse_version('2.5.0'):
+        crs = pyproj.CRS.from_epsg(epsg_number)
+    else:
+        import fiona
+        crs = fiona.crs.from_epsg(epsg_number)
+    return crs
+
+
 def getCRSfromRasterFile(raster_file):
     """Create a CRS object from a raster file.
 
@@ -254,22 +283,20 @@ def getUTMfromWGS84(longitude, latitude):
     def get_utm_zone(longitude):
         return int(1 + (longitude + 180.0) / 6.0)
 
-    def is_northern(latitude):
-        """
-        Determines if given latitude is a northern for UTM
-        """
-        if latitude < 0.0:
-            return 0
-        else:
-            return 1
-
-    utm_crs = osr.SpatialReference()
-    utm_crs.SetWellKnownGeogCS("WGS84")  # Set geographic coordinate system to handle latitude/longitude
     zone = get_utm_zone(longitude)
-    utm_crs.SetUTM(get_utm_zone(longitude), is_northern(latitude))
-    wgs84_crs = utm_crs.CloneGeogCS()  # Clone ONLY the geographic coordinate system
 
-    return zone, utm_crs, wgs84_crs
+    # Determines if given latitude is a northern for UTM
+    is_northern = int(latitude > 0.0)
+
+    osr_srs = osr.SpatialReference()
+    osr_srs.SetWellKnownGeogCS("WGS84")  # Set geographic coordinate system to handle latitude/longitude
+    osr_srs.SetUTM(zone, is_northern)
+    osr_srs.AutoIdentifyEPSG()
+    out_epsg = int(osr_srs.GetAuthorityCode(None))
+    osr_srs.ImportFromEPSG(out_epsg)
+    wgs84_crs = osr_srs.CloneGeogCS()  # Clone ONLY the geographic coordinate system
+
+    return zone, osr_srs, wgs84_crs, out_epsg
 
 
 def getProjectedCRSForXY(x_coord, y_coord, xy_epsg=4326):
@@ -298,46 +325,41 @@ def getProjectedCRSForXY(x_coord, y_coord, xy_epsg=4326):
 
         """
 
-    # Based On :https://stackoverflow.com/a/10239676
-    from pyproj import Proj, transform
-
     # Coordinates need to be in wgs84 so project them
     if xy_epsg != 4326:
-        inProj = Proj(init='epsg:{}'.format(xy_epsg))
-        outProj = Proj(init='epsg:4326')
-        longitude, latitude = transform(inProj, outProj, x_coord, y_coord)
-    else:
-        longitude, latitude = x_coord, y_coord
+        if parse_version(pyproj.__version__) >= parse_version('2.5.0'):
+            inProj = from_epsg(xy_epsg)                  #Proj(init='epsg:{}'.format(xy_epsg))
+            outProj = from_epsg(4326)                    #Proj(init='epsg:4326')
 
-    utm_zone = int(1 + (longitude + 180.0) / 6.0)
+            from pyproj import Transformer
+            transformer = Transformer.from_crs(inProj, outProj,always_xy=True)
+            longitude, latitude = transformer.transform(x_coord, y_coord)
+        else:
+            # Based On :https://stackoverflow.com/a/10239676
+            from pyproj import Proj, transform
+            inProj =  Proj(init='epsg:{}'.format(xy_epsg))
+            outProj = Proj(init='epsg:4326')
+            longitude, latitude = transform(inProj, outProj, x_coord, y_coord)
 
-    # Determines if given latitude is a northern for UTM        1 is northern, 0 is southern
-    is_northern = int(latitude > 0.0)
+    else:
+        longitude, latitude = x_coord, y_coord
 
-    utm_crs = crs()
-    utm_crs.srs = osr.SpatialReference()
+    utm_zone, _, _, out_epsg = getUTMfromWGS84(longitude, latitude)
 
+    osr_srs = osr.SpatialReference()
     # Use Australian or New Zealand local systems otherwise use global wgs84 UTM zones.
     if (108.0 <= longitude <= 155.0) and (-45.0 <= latitude <= -10.0):
         # set to Australian GDA94 MGA Zone XX
-        utm_crs.srs.ImportFromEPSG(int('283{}'.format(utm_zone)))
+        out_epsg = int('283{}'.format(utm_zone))
     elif (166.33 <= longitude <= 178.6) and (-47.4 <= latitude <= -34.0):
         # set to NZGD2000 / New Zealand Transverse Mercator 2000
-        utm_crs.srs.ImportFromEPSG(2193)
-    else:
-        # Set to Global WGS84 Utm Zone.
-        utm_crs.srs.SetWellKnownGeogCS("WGS84")
-        utm_crs.srs.SetUTM(utm_zone, is_northern)
-        utm_crs.srs.AutoIdentifyEPSG()
-
-    utm_crs.set_epsg(utm_crs.srs.GetAuthorityCode(None))
+        out_epsg=2193
 
-    utm_crs.epsg_predicted = False
-    utm_crs.crs_wkt = utm_crs.srs.ExportToWkt()
-    utm_crs.proj4 = utm_crs.srs.ExportToProj4()
+    utmzone_crs = crs()
+    utmzone_crs.getFromEPSG(out_epsg)
 
-    if utm_crs.srs.IsProjected():
-        return utm_crs
+    if utmzone_crs.srs.IsProjected():
+        return utmzone_crs
     else:
         return None
 
@@ -364,7 +386,7 @@ def distance_metres_to_dd(longitude, latitude, distance_metres):
         if not isinstance(argCheck[1], float):
             raise TypeError('{} must be a floating number.'.format(argCheck[0]))
 
-    if not isinstance(distance_metres, (int, long, float)):
+    if not isinstance(distance_metres, six.integer_types + (float, )):
         raise TypeError('distance_metres must be a floating number.')
 
     # get the required Spatial reference systems
@@ -374,14 +396,22 @@ def distance_metres_to_dd(longitude, latitude, distance_metres):
 
     # create transform component
     wgs842utm_transform = osr.CoordinateTransformation(wgs84SRS, utm_crs.srs)  # (<from>, <to>)
-    # Project to UTM
-    easting, northing, alt = wgs842utm_transform.TransformPoint(longitude, latitude, 0)
 
     # create transform component
     utm2wgs84_transform = osr.CoordinateTransformation(utm_crs.srs, wgs84SRS)  # (<from>, <to>)
 
-    # Project back to WGS84 after adding distance to eastings. returns easting, northing, altitude
-    newLong, newLat, alt = utm2wgs84_transform.TransformPoint(easting + distance_metres, northing, 0)
+    # Project to UTM
+    easting, northing, alt = wgs842utm_transform.TransformPoint(longitude, latitude, 0)
+
+    # if easting/northing return infity, then swap lat and long over.
+    if easting == float("inf"):
+        easting, northing, alt = wgs842utm_transform.TransformPoint(latitude,longitude, 0)
+
+        # Project back to WGS84 after adding distance to eastings. returns easting, northing, altitude
+        newLat, newLong, alt = utm2wgs84_transform.TransformPoint(easting + distance_metres, northing, 0)
+    else:
+        # Project back to WGS84 after adding distance to eastings. returns easting, northing, altitude
+        newLong, newLat, alt = utm2wgs84_transform.TransformPoint(easting + distance_metres, northing, 0)
 
     # Create a point objects.
     origPoint = geometry.Point(longitude, latitude)
@@ -395,7 +425,3 @@ def distance_metres_to_dd(longitude, latitude, distance_metres):
         distance_dd = -distance_dd
 
     return distance_dd
-
-
-
-
diff --git a/pyprecag/describe.py b/pyprecag/describe.py
index 36615b4..31e3775 100644
--- a/pyprecag/describe.py
+++ b/pyprecag/describe.py
@@ -1,6 +1,7 @@
 import csv
 import datetime
 import difflib
+import io
 import logging
 import os
 import re
@@ -14,6 +15,7 @@
 import geopandas
 import numpy as np
 import pandas as pd
+import six
 from geopandas import GeoDataFrame
 from osgeo import gdal
 
@@ -56,6 +58,12 @@ def __init__(self, input_data):
 
         return
 
+    @staticmethod
+    def get_character_encoding(filename):
+        with io.open(filename, 'rb') as f:
+           data = f.read(1024 * 100)
+        return chardet.detect(data)['encoding']
+
     def open_geo_dataframe(self):
         """Create geopandas from file"""
         return GeoDataFrame.from_file(self.source, encoding=self.file_encoding)
@@ -65,17 +73,13 @@ def describeFile(self):
         Describe a vector File and set class properties
         """
 
+        self.file_encoding = self.get_character_encoding(self.source)
         # Use this open so as to not hold open and use up memory
         gdf = GeoDataFrame.from_file(self.source, encoding=self.file_encoding)
 
         with fiona.open(self.source) as fio_coll:
             self.crs.getFromWKT(fio_coll.crs_wkt)
 
-        # No idea why this works but it does so use it.
-        exec ('rawstring = "{}"'.format(repr(','.join(gdf.columns.values))))
-        result = chardet.detect(rawstring)
-        self.file_encoding = result['encoding']
-
         self.feature_count = len(gdf)
         self.extent = list(gdf.total_bounds)
 
@@ -147,7 +151,7 @@ def open_pandas_dataframe(self, **kwargs):
         return pdf
 
     def get_column_names(self):
-        return self.column_properties.keys()
+        return list(self.column_properties.keys())
 
     def get_alias_column_names(self):
         return [val['alias'] for key, val in self.column_properties.items()]
@@ -158,32 +162,28 @@ def get_shapefile_column_names(self):
     def get_column_types(self):
         return [val['type'] for key, val in self.column_properties.items()]
 
+    @staticmethod
+    def get_character_encoding(filename):
+        with io.open(filename, 'rb') as f:
+           data = f.read(1024 * 100)
+        return chardet.detect(data)['encoding']
+
     def describe_file(self):
         """Describe a CSV File and set class properties
         """
-        with open(self.source, 'r') as f:
-            # sniff into 10KB of the file to check its dialect
-            # this will sort out the delimiter and quote character.
-            self.dialect = csv.Sniffer().sniff(f.read(10 * 1024))
-            f.seek(0)  # reset read to start of file
-
-            # read header based on the 10k of file.
-            header = csv.Sniffer().has_header(f.read(10 * 1024))
-            f.seek(0)  # reset read to start of file
-            if not header:
-                warnings.warn("The CSV file doesn't appear to contain column headers")
-                self.has_column_header = False
-
-            f.seek(0)  # reset read to start of file
-
-        detector = chardet.UniversalDetector()
-        with open(self.source, 'rb') as eaop:
-            for line in eaop.readlines(100):
-                detector.feed(line)
-                if detector.done:
-                    break
-            detector.close()
-        self.file_encoding = detector.result['encoding']
+
+        self.file_encoding = self.get_character_encoding(self.source)
+
+        # reopen the file as text, using the detected encoding
+        with io.open(self.source, mode='rt', encoding=self.file_encoding) as f:
+            text = six.ensure_str(f.read(10240))
+        self.dialect = csv.Sniffer().sniff(text)
+        # read header based on the decoded text
+        header = csv.Sniffer().has_header(text)
+        if not header:
+            warnings.warn("The CSV file doesn't appear to contain column headers")
+            self.has_column_header = False
+
 
         pandas_df = self.open_pandas_dataframe()
         self.row_count = len(pandas_df)
@@ -282,12 +282,6 @@ def save_geopandas_tofile(inputGeoDataFrame, output_filename, overwrite=True, fi
                                               dur=datetime.timedelta(seconds=time.time() - step_time)))
 
 
-def get_dataframe_encoding(dataframe):
-    exec ('rawstring = "{}"'.format(repr(','.join(dataframe.columns))))
-    result = chardet.detect(rawstring)
-    return result['encoding']
-
-
 def get_column_properties(dataframe):
     """ Get a dictionary representing Column Properties for a pandas dataframe or a geopandas geodataframe.
        Includes:
@@ -321,14 +315,18 @@ def get_column_properties(dataframe):
             if fldtype == 'long':
                 fldtype = 'int'
 
-        if isinstance(col, unicode):
-            aliasFld = unidecode(unicode(col))
-        else:
-            aliasFld = col
+        aliasFld = col
+        # check for unicode characters
+        if not all(ord(char) < 128 for char in col):
+            try: # python 2.7
+                if isinstance( col, unicode):
+                    aliasFld = unidecode(col)
+            except: # python 3.7
+                aliasFld = unidecode(six.ensure_str(col))
 
         # create a shapefile valid name 10 alpha numeric and underscore characters.
         # to keep underscore('_'), addit after the 9
-        shpFld = re.sub('[^A-Za-z0-9_-]+', '', col)[:10]
+        shpFld = re.sub('[^A-Za-z0-9_-]+', '', aliasFld)[:10]
 
         column_desc[col] = {'alias': aliasFld.replace(' ', ''),
                             'shapefile': shpFld,
@@ -366,12 +364,16 @@ def predictCoordinateColumnNames(column_names):
                     seqMatchDict[guess] = difflib.SequenceMatcher(None, eaFld.upper(), guess.upper(), True).ratio()
 
                 # create short list of matches and ratios
-                valList.append(max(seqMatchDict.iteritems(), key=lambda x: x[1]))
+                iter_items = six.iteritems(seqMatchDict)
+                valList.append(max(iter_items, key=lambda x: x[1]))
 
         # select the largest ratio as the best match
         if len(valList) > 0:
             best_match = max(valList, key=itemgetter(1))[0]
-            exec ('{}_column = "{}"'.format(eaVal, best_match))
+            if eaVal == 'x':
+                x_column = best_match
+            else:
+                y_column = best_match
 
     LOGGER.debug('GeoCSV Columns:     x = {}, y = {}'.format(x_column, y_column))
     return [x_column, y_column]
diff --git a/pyprecag/kriging_ops.py b/pyprecag/kriging_ops.py
index 3fa3636..11ccc9d 100644
--- a/pyprecag/kriging_ops.py
+++ b/pyprecag/kriging_ops.py
@@ -1,5 +1,5 @@
 import collections
-import datetime
+from datetime import datetime, timedelta
 import glob
 import inspect
 import logging
@@ -15,14 +15,15 @@
 import pandas as pd
 import rasterio
 from rasterio import features
+import six
 
-from unidecode import unidecode
 from collections import OrderedDict
 
 from . import config
 from .convert import add_point_geometry_to_dataframe, numeric_pixelsize_to_string
 from .describe import predictCoordinateColumnNames
 from .raster_ops import raster_snap_extent, create_raster_transform
+from pyprecag import number_types
 
 LOGGER = logging.getLogger(__name__)
 LOGGER.addHandler(logging.NullHandler())
@@ -39,17 +40,18 @@ def test_for_windows():
 
 # consider enums - https://stackoverflow.com/questions/36932/how-can-i-represent-an-enum-in-python
 # The codes for the variogram. The key is the tag used in the control file_csv
-VESPER_OPTIONS = {"jigraph": {"Don't Show": 0, "Show": 1},
-                  "jimap": {"Don't Show": 0, "Show": 1},
-                  "jlockrg": {"Local": 1, "Global": 0},
-                  "jpntkrg": {"Block": 0, "Point": 1, "Punctual": 1},
-                  "jsetrad": {"Calculate Radius": 0, "Set Radius": 1},
-                  "jcomvar": {"Define Variogram Parameter": 0, "Compute Variogram": 1},
-                  "modtyp": {"Spherical": 1, "Exponential": 2, "Gaussian": 3,
-                             "Linear with sill": 4, "Stable": 5, "Generalised Cauchy": 6,
-                             "Matern": 7, "Double spherical": 8, "Double exponential": 9},
-                  "iwei": {"unity": 0, "No. of pairs": 1, "1/variance": 2, "no_pairs/variance": 3,
-                           "no_pairs/fitted": 4}}
+VESPER_OPTIONS = {"jigraph": OrderedDict([("Don't Show", 0), ("Show", 1)]),
+                  "jimap": OrderedDict([("Don't Show", 0), ("Show", 1)]),
+                  "jlockrg": OrderedDict([("Local", 1), ("Global", 0)]),
+                  "jpntkrg": OrderedDict([("Punctual", 1), ("Block", 0), ("Point", 1)]),
+                  "jsetrad": OrderedDict([("Calculate Radius", 0), ("Set Radius", 1)]),
+                  "jcomvar": OrderedDict([("Define Variogram Parameter", 0), ("Compute Variogram", 1)]),
+                  "modtyp": OrderedDict([("Spherical", 1), ("Exponential", 2), ("Gaussian", 3),
+                                         ("Linear with sill", 4), ("Stable", 5), ("Generalised Cauchy", 6),
+                                         ("Matern", 7), ("Double spherical", 8), ("Double exponential", 9) ]),
+                  "iwei": OrderedDict([("unity", 0), ("No. of pairs", 1), ("1/variance", 2),
+                                       ("no_pairs/variance", 3), ("no_pairs/fitted", 4)])
+                  }
 
 
 class VesperControl(collections.MutableMapping, dict):
@@ -102,20 +104,20 @@ def __setitem__(self, key, value):
             raise AttributeError('VesperControl has no attribute {}'.format(key))
 
         if key in VESPER_OPTIONS.keys():
-            if isinstance(value, basestring):
+            if isinstance(value, six.string_types):
                 try:
                     value = VESPER_OPTIONS[key][value]
                 except KeyError:
                     raise ValueError('{} is an invalid option for {}. Options are {}'.format(
-                        value, key, VESPER_OPTIONS[key]))
+                        value, key, dict(VESPER_OPTIONS[key])))
 
-            elif isinstance(value, (int, float, long)):
+            elif isinstance(value, number_types):
                 if value not in VESPER_OPTIONS[key].values():
                     raise ValueError('{} is an invalid option for {}. Options are {}'.format(
-                        value, key, VESPER_OPTIONS[key]))
+                        value, key, dict(VESPER_OPTIONS[key])))
 
         if isinstance(self.__defaults[key], float):
-            if not isinstance(value, (int, long, float)):
+            if not isinstance(value, number_types):
                 raise ValueError('value for key {} should be a float or int - Got {}'.format(
                     key, type(value.__name__)))
 
@@ -192,7 +194,7 @@ def vesper_text_to_raster(control_textfile, krig_epsg=0, nodata_value=-9999):
         raise IOError("Invalid path: {}".format(control_textfile))
 
     for argCheck in [('krig_epsg', krig_epsg), ('nodata_value', nodata_value)]:
-        if not isinstance(argCheck[1], (int, long)):
+        if not isinstance(argCheck[1], six.integer_types):
             raise TypeError('{} must be a integer.'.format(argCheck[0]))
 
     start_time = time.time()
@@ -219,7 +221,7 @@ def vesper_text_to_raster(control_textfile, krig_epsg=0, nodata_value=-9999):
         ci_file.writelines("Median Prediction SE    : {:.5f}\n".format(median_val))
         ci_file.writelines("95% Confidence Interval : {:.5f}\n\n".format(2 * 1.96 * median_val))
         ci_file.writelines(
-            "Date/time : " + datetime.datetime.now().strftime("%d-%b-%Y %H:%M") + "\n")
+            "Date/time : " + datetime.now().strftime("%d-%b-%Y %H:%M") + "\n")
         userhome = os.path.expanduser('~')
         ci_file.writelines("Username  : " + os.path.split(userhome)[-1] + "\n")
 
@@ -271,16 +273,15 @@ def vesper_text_to_raster(control_textfile, krig_epsg=0, nodata_value=-9999):
                                     transform=out_se.transform, fill=nodata_value)
         out_se.write(burned, indexes=1)
 
-    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                               dur=datetime.timedelta(
-                                                   seconds=time.time() - start_time)))
+    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(
+        inspect.currentframe().f_code.co_name, '',
+        dur=str(timedelta( seconds=time.time() - start_time))))
     return out_pred_tif, out_se_tif, out_ci_txt
 
 
 def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folder,
                              control_textfile='', coord_columns=[], epsg=0, display_graphics=False,
                              control_options=VesperControl(),
-                             block_size=10,
                              vesper_exe=vesper_exe):
     """Prepare data for vesper kriging and create a windows batch file to run outside the
     python/pyprecag environment.
@@ -290,7 +291,7 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
     A description of the keys and their options can be found in the VESPER user manual at
     https://sydney.edu.au/agriculture/pal/documents/Vesper_1.6_User_Manual.pdf
 
-    block_size will be replaced by the xside and yside keys in the VesperControl Object.
+    block_size has been removed from version 1.0.0.  Please use use VesperControl xsize and ysize instead.
 
     Outputs:  The following files will be added to the vesper sub-folder in the output folder.
         *_control_*.txt  - The vesper control file
@@ -311,7 +312,7 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
         grid_filename (str): The vesper grid file.
         control_textfile (str): The name of the control text file without the path
         out_folder (str): The folder to add outputs too. A 'Vesper' sub directory will be created
-        block_size (int): The size to apply for block Kriging. Units are from the coordinate system.
+
         coord_columns (List): The columns representing the X and Y coordinates.
         epsg (int) : The epsg_number number for the data. If 0 the en_epsg or
                      enepsg column (if exists) will be used.
@@ -322,21 +323,12 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
        vesper_batfile, vesper_ctrlfile: The paths to the generated batch file and control file.
     """
 
-    # Vesper only works with Windows
-    # test_for_windows()
-
     if not isinstance(in_dataframe, (gpd.GeoDataFrame, pd.DataFrame)):
         raise TypeError('Invalid input data :in_dataframe')
 
     if not os.path.exists(grid_filename):
         raise IOError("Invalid path: {}".format(grid_filename))
 
-    if not isinstance(block_size, (int, long)):
-        raise TypeError('block_size must be an integer'.format(block_size))
-
-    warnings.warn('block_size is deprecated, use VesperControl xsize and ysize instead',
-                  PendingDeprecationWarning)
-
     if out_folder.strip() in [None, '']:
         raise TypeError('Please specify an output folder')
 
@@ -353,7 +345,7 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
         if eaFld not in in_dataframe.columns:
             raise TypeError('Column {} does not exist'.format(eaFld))
 
-    if not isinstance(epsg, (int, long)):
+    if not isinstance(epsg, six.integer_types):
         raise TypeError('EPSG {} must be a integer.'.format(epsg))
 
     if not isinstance(control_options, VesperControl):
@@ -362,7 +354,7 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
     start_time = time.time()
 
     # Create a filename compatible copy of the krig_column
-    krig_col_file = re.sub('[^A-Za-z0-9_-]+', '', unidecode(unicode(krig_column)))
+    krig_col_file = re.sub('[^A-Za-z0-9_-]+', '', six.ensure_str(krig_column, encoding='ascii', errors='ignore'))
 
     out_sub_name = os.path.basename(grid_filename)[:20]
 
@@ -403,8 +395,8 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
     # in use either by vesper or another app.
 
     # Start always start with the control file. this is what gets used within QGIS
-    files_list = glob.glob(
-        os.path.join(vesper_outdir, "{}_*_{}.*".format(out_sub_name, krig_col_file)))
+    files_list = glob.glob(os.path.join(vesper_outdir, "{}_*_{}.*".format(out_sub_name, krig_col_file)))
+
     if vesper_ctrlfile in files_list:
         files_list.insert(0, files_list.pop(files_list.index(vesper_ctrlfile)))
 
@@ -417,8 +409,7 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
 
     # if control file already exists then replace it, and delete all matching kriging
     # outputs and tiff files.
-    for eaFile in glob.glob(
-            os.path.join(vesper_outdir, "{}_*_{}.*".format(out_sub_name, krig_col_file))):
+    for eaFile in glob.glob(os.path.join(vesper_outdir, "{}_*_{}.*".format(out_sub_name, krig_col_file))):
         os.remove(eaFile)
         LOGGER.debug('Deleted file {}'.format(eaFile))
 
@@ -461,11 +452,6 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
                            jimap=int(display_graphics),  # 1, show map, otherwise 0)
                            )
 
-    # high density kriging. fix for pre VesperControl Class
-    if block_size != 10:  # only update if not the default variable value
-        control_options.update({"xside": block_size,
-                                "yside": block_size})
-
     # ---------------------------------------------------------------------------------------------
     # Write a VESPER control file
     control_options.write_to_file(vesper_ctrlfile)
@@ -505,9 +491,10 @@ def prepare_for_vesper_krige(in_dataframe, krig_column, grid_filename, out_folde
         with open(vesper_batfile, 'w') as wBatFile:
             wBatFile.write(bat_file_string)
 
-    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                               dur=datetime.timedelta(
-                                                   seconds=time.time() - start_time)))
+    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(
+        inspect.currentframe().f_code.co_name, '',
+        dur=str(timedelta(seconds=time.time() - start_time))
+    ))
     return vesper_batfile, vesper_ctrlfile
 
 
@@ -542,5 +529,4 @@ def run_vesper(ctrl_file, bMinimiseWindow=True, vesper_exe=vesper_exe):
     process = subprocess.Popen([vesper_exe, ctrl_file], cwd=os.path.dirname(ctrl_file),
                                startupinfo=info)
     process.wait()
-    logging.info('{:<30}\t{dur:<15}'.format('Vesper Kriging', dur=datetime.timedelta(
-        seconds=time.time() - task_time)))
+    LOGGER.info('{:<30}\t{dur:<15}'.format('Vesper Kriging', dur=str(timedelta(seconds=time.time() - task_time))))
diff --git a/pyprecag/processing.py b/pyprecag/processing.py
index 14eea94..87f2266 100644
--- a/pyprecag/processing.py
+++ b/pyprecag/processing.py
@@ -4,9 +4,11 @@
 import re
 from collections import defaultdict
 
-from itertools import izip
 import logging
 import os
+import six
+from rasterio import crs as riocrs
+from six.moves import zip as izip
 
 import random
 import time
@@ -18,7 +20,8 @@
 
 import rasterio
 
-from fiona.crs import from_epsg
+from pyprecag.crs import from_epsg
+
 from geopandas import GeoDataFrame, GeoSeries
 from osgeo import gdal
 
@@ -27,7 +30,7 @@
 from rasterio.fill import fillnodata
 from rasterio.mask import mask as rio_mask
 from rasterio.warp import reproject, Resampling
-# from rasterio.warp import aligned_target
+
 from rasterio.windows import get_data_window, intersection, from_bounds
 from scipy.cluster.vq import *
 from scipy import stats
@@ -43,69 +46,123 @@
     convert_polygon_feature_to_raster, drop_z, deg_to_8_compass_pts, point_to_point_bearing, \
     text_rotation
 
-from .describe import save_geopandas_tofile, VectorDescribe, get_dataframe_encoding
+from .describe import save_geopandas_tofile, VectorDescribe
 from .errors import GeometryError, SpatialReferenceError
 from .vector_ops import thin_point_by_distance
 from .raster_ops import focal_statistics, save_in_memory_raster_to_file, reproject_image, \
     calculate_image_indices, stack_and_clip_rasters
 
+from pyprecag import number_types
+
 LOGGER = logging.getLogger(__name__)
 LOGGER.addHandler(logging.NullHandler())
 
 
-def block_grid(in_shapefilename, pixel_size, out_rasterfilename,
-               out_vesperfilename, nodata_val=-9999, snap=True, overwrite=False):
+def block_grid(in_shapefilename, pixel_size, out_rasterfilename, out_vesperfilename,
+              out_epsg=0,groupby=None, nodata_val=-9999, snap=True, overwrite=False):
     """Convert a polygon boundary to a 0,1 raster and generate a VESPER compatible list of
        coordinates for kriging.
 
-        Args:
+       The input polygon shapefile will be reproject to the out_epsg coordinate system.
 
+       Will process in a loop based on the groupby field name. All polygons with the same
+       attribute will be dissolved together and a single block grid will be created. The
+       groupby attribute will be added to the output filename.
+
+        Args:
             in_shapefilename (str): Input polygon shapefile
             pixel_size (float):  The required output pixel size
             out_rasterfilename (str): Filename of the raster Tiff that will be created
             out_vesperfilename (str):  The output vesper file
+            out_epsg (int): The epsg number for the output raster coordinate system.
+            groupby (str): Column name to group polygon features with.
             nodata_val (int): an integer to use as nodata
             snap (bool): Snap Extent to a factor of the Pixel size
             overwrite (bool): if true overwrite existing file
 
-        Requirements: Input shapefile should be in a projected coordinate system........
-
-        Notes:
-            Define pixel_size and no data value of new raster
-            see http://stackoverflow.com/questions/2220749/rasterizing-a-gdal-layer
-                https://gis.stackexchange.com/a/31658
-
-            This works, but the extent of the output differs from that of an arcpy generated raster.
-            See post: https://gis.stackexchange.com/q/139336
     """
 
-    if not isinstance(pixel_size, (int, long, float)):
+    if not isinstance(pixel_size, number_types):
         raise TypeError('Pixel size must be an integer or floating number.')
 
-    if not isinstance(nodata_val, (int, long)):
-        raise TypeError('Nodata value must be an integer.')
+    
+    for ea_arg in [('nodata_val', nodata_val), ('out_epsg', out_epsg)]:
+        if not isinstance(ea_arg[1], six.integer_types):
+            raise TypeError('{} must be a integer - Got {}'.format(*ea_arg))
+        
+    if out_epsg <= 0:
+        raise ValueError('EPSG Code ({}) must be a positive integer above 0'.format(out_epsg))
+
+    if riocrs.CRS.from_epsg(out_epsg).is_geographic:
+        raise ValueError('EPSG Code ({}) must be a projected coordinate system'.format(out_epsg))
 
     desc_poly_shp = VectorDescribe(in_shapefilename)
 
-    if not desc_poly_shp.crs.srs.IsProjected():
-        raise SpatialReferenceError('Shapefile must be in a projected coordinate system')
+    gdf_poly = desc_poly_shp.open_geo_dataframe()
 
     if 'POLY' not in desc_poly_shp.geometry_type.upper():
         raise GeometryError('Invalid Geometry. Input shapefile should be polygon or multipolygon')
 
-    start_time = time.time()
-    convert_polygon_to_grid(in_shapefilename=in_shapefilename,
-                            out_rasterfilename=out_rasterfilename,
-                            pixel_size=pixel_size,
-                            snap_extent_to_pixel=snap,
-                            nodata_val=nodata_val,
-                            overwrite=overwrite)
+    if groupby is None or groupby.strip() == '':
+        groupby = None
 
-    convert_grid_to_vesper(in_rasterfilename=out_rasterfilename,
-                           out_vesperfilename=out_vesperfilename)
+    if groupby is not None and groupby not in gdf_poly.columns:
+        raise ValueError('Groupby column {} does not exist'.format(groupby))
 
-    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                               dur=timedelta(seconds=time.time() - start_time)))
+    #-------------------------------------------------------------------------------------------
+    # reproject shapefile
+    if out_epsg != 0 and  desc_poly_shp.crs.epsg_number != out_epsg:
+        gdf_poly.to_crs(epsg=out_epsg,inplace=True)
+
+    del desc_poly_shp
+
+    # If a column name is specified, dissolve by this and create multi-polygons
+    # otherwise dissolve without a column name at treat all features as one
+    if groupby is not None:
+        # change null/nones to blank string
+        gdf_poly[groupby].fillna('', inplace=True)
+
+        gdf_poly = gdf_poly.dissolve(groupby, as_index=False).copy()
+
+    else:
+        gdf_poly = GeoDataFrame(geometry=[gdf_poly.unary_union])
+
+    # Loop through polygons features ----------------------------------------------------------------
+    output_files = []
+    for index, feat in gdf_poly.iterrows():
+        loop_time = time.time()
+        step_time = time.time()
+        status = '{} of {}'.format(index + 1, len(gdf_poly))
+
+        if groupby is not None:
+            feat_name = re.sub('[^0-9a-zA-Z]+', '-', str(feat[groupby])).strip('-')
+            if feat_name == '':
+                feat_name = 'No-Name'
+
+            r_file, r_ext = os.path.splitext(out_rasterfilename)
+            out_file = r_file + "_" + feat_name + r_ext
+        else:
+            out_file = out_rasterfilename
+            feat_name = 'All Polygons'
+
+        # From shapes create a blockgrid mask. -----------------------------------------------------
+        blockgrid, new_blockmeta = convert_polygon_feature_to_raster(feat, pixel_size)
+
+        with rasterio.open(out_file, 'w', driver='GTiff', count=1,tfw='YES',
+                           crs=from_epsg(out_epsg), **new_blockmeta) as dest:
+            dest.write(blockgrid, 1)
+        output_files += [out_file]
+
+        del blockgrid, new_blockmeta
+
+        LOGGER.info('{:<30} {:>10}   {:<15} {dur}'.format('Feature to block_grid', status, feat_name,
+                                                          dur=timedelta(seconds=time.time() - step_time)))
+
+        r_file, r_ext = os.path.splitext(out_file)
+        convert_grid_to_vesper(in_rasterfilename=out_file,
+                           out_vesperfilename=out_file.replace(r_ext, '_v.txt'))
+
+    return output_files
 
 
 def create_polygon_from_point_trail(points_geodataframe, points_crs, out_filename, thin_dist_m=1.0,
@@ -148,7 +205,7 @@ def create_polygon_from_point_trail(points_geodataframe, points_crs, out_filenam
 
     for argCheck in [('thin_dist_m', thin_dist_m), ('aggregate_dist_m', aggregate_dist_m),
                      ('buffer_dist_m', buffer_dist_m), ('shrink_dist_m', shrink_dist_m)]:
-        if not isinstance(argCheck[1], (int, long, float)):
+        if not isinstance(argCheck[1], number_types):
             raise TypeError('{} must be a floating number.'.format(argCheck[0]))
 
     if not isinstance(points_geodataframe, GeoDataFrame):
@@ -275,8 +332,10 @@ def create_polygon_from_point_trail(points_geodataframe, points_crs, out_filenam
 
     save_geopandas_tofile(gdf_final, out_filename, overwrite=True)
 
-    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                               dur=timedelta(seconds=time.time() - start_time)))
+    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(
+        inspect.currentframe().f_code.co_name, '',
+        dur=str(timedelta(seconds=time.time() - start_time))
+    ))
 
     thin_ratio = (4 * 3.14 * gdf_final['Area'].sum() /
                   (gdf_final['Perimeter'].sum() * gdf_final['Perimeter'].sum()))
@@ -352,7 +411,7 @@ def clean_trim_points(points_geodataframe, points_crs, process_column, output_cs
             raise IOError('Output folder does not exist: {}'.format(os.path.dirname(argCheck)))
 
     for argCheck in [('thinDist_m', thin_dist_m), ('stdev', stdevs)]:
-        if not isinstance(argCheck[1], (int, long, float)):
+        if not isinstance(argCheck[1], number_types):
             raise TypeError('{} must be a integer or floating number.'.format(argCheck[0]))
 
     if process_column not in points_geodataframe.columns:
@@ -534,6 +593,7 @@ def clean_trim_points(points_geodataframe, points_crs, process_column, output_cs
     gdf_points['Northing'] = gdf_points.geometry.apply(lambda p: p.y)
     gdf_points['EN_EPSG'] = points_crs.epsg_number
     gdf_points.crs = points_crs.epsg
+    
     # Clean up the original input dataframe and remove existing geometry and coord columns
     alt_coord_columns = config.get_config_key('geoCSV')['xCoordinate_ColumnName']
     alt_coord_columns += config.get_config_key('geoCSV')['yCoordinate_ColumnName']
@@ -552,7 +612,12 @@ def clean_trim_points(points_geodataframe, points_crs, process_column, output_cs
     gdf_final.drop([id_col], inplace=True, axis=1)
 
     # get the appropriate file encoding and save the file
-    file_encoding = get_dataframe_encoding(gdf_final)
+    # TODO: CHECK THIS I dont think setting a default in the config will work.
+    try:
+        file_encoding = config.read_config()['geoCSV']['file_encoding']
+    except KeyError:
+        file_encoding = 'utf-8'
+    #file_encoding = get_dataframe_encoding(gdf_final)
 
     gdf_final[gdf_final['filter'].isnull()].drop(['geometry', 'filter'], axis=1) \
         .to_csv(output_csvfile, index=False, encoding=file_encoding)
@@ -587,8 +652,10 @@ def clean_trim_points(points_geodataframe, points_crs, process_column, output_cs
     LOGGER.info('\nResults:---------------------------------------\n{}\n'.format(
         results_table.to_string(index=False, justify='center')))
 
-    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(inspect.currentframe().f_code.co_name, '',
-                                               dur=timedelta(seconds=time.time() - start_time)))
+    LOGGER.info('{:<30}\t{dur:<15}\t{}'.format(
+        inspect.currentframe().f_code.co_name, '',
+        dur=str(timedelta(seconds=time.time() - start_time)))
+    )
 
     return gdf_final[gdf_final['filter'].isnull()], points_crs
 
@@ -619,7 +686,7 @@ def random_pixel_selection(raster, raster_crs, num_points, out_shapefile=None):
         raise TypeError("Input should be a rasterio.DatasetReader created "
                         "using rasterio.open(os.path.normpath())")
 
-    if not isinstance(num_points, (int, long)):
+    if not isinstance(num_points, six.integer_types):
         raise TypeError('Size must be an Integer.')
 
     if not isinstance(raster_crs, pyprecag_crs.crs):
@@ -844,12 +911,14 @@ def extract_pixel_statistics_for_points(points_geodataframe, points_crs, rasterf
         points_geodataframe.drop(['geometry'], axis=1).to_csv(output_csvfile, index=False)
 
         LOGGER.info('{:<30}\t{:>10}   {dur:<15} {}'.format(
-            'Saved CSV', '', os.path.basename(output_csvfile),
-            dur=timedelta(seconds=time.time() - step_time)))
+                'Saved CSV', '', os.path.basename(output_csvfile),
+                dur=str(timedelta(seconds=time.time() - step_time))
+        ))
 
     LOGGER.info('{:<30}\t{:>10} {dur:>15}\t{}'.format(
         inspect.currentframe().f_code.co_name, '', '',
-        dur=timedelta(seconds=time.time() - start_time)))
+        dur=str(timedelta(seconds=time.time() - start_time))
+        ))
 
     return points_geodataframe, points_crs
 
@@ -916,10 +985,10 @@ def multi_block_bands_processing(image_file, pixel_size, out_folder, band_nums=[
             raise IOError('{} does not exist-Got {}'.format(ea_arg[0], os.path.dirname(ea_arg[1])))
 
     for ea_arg in [('pixel_size', pixel_size), ('image_nodata', image_nodata)]:
-        if not isinstance(ea_arg[1], (int, long, float)):
+        if not isinstance(ea_arg[1], number_types):
             raise TypeError('{} must be a integer or floating number - Got {}'.format(*ea_arg))
 
-    if not isinstance(image_epsg, (int, long)):
+    if not isinstance(image_epsg, six.integer_types):
         raise TypeError('image_epsg must be a integer - Got {}'.format(image_epsg))
 
     if not isinstance(band_nums, list):
@@ -1005,7 +1074,7 @@ def multi_block_bands_processing(image_file, pixel_size, out_folder, band_nums=[
         poly_epsg = desc_poly.crs.epsg_number
         if poly_epsg != image_epsg:
             # we have to reproject the vector points to match the raster
-            print 'WARNING: Projecting points from {} to {}'.format(desc_poly.crs.epsg, image_epsg)
+            print('WARNING: Projecting points from {} to {}'.format(desc_poly.crs.epsg, image_epsg))
             gdf_poly.to_crs(epsg=image_epsg, inplace=True)
             poly_epsg = image_epsg
 
@@ -1078,7 +1147,7 @@ def multi_block_bands_processing(image_file, pixel_size, out_folder, band_nums=[
                 for iband in range(1, src.count + 1):
                     # image statistics have changed so don't copy the tags
                     cleaned_tags = dict(
-                        [(key, val) for key, val in src.tags(iband).iteritems() if
+                        [(key, val) for key, val in six.iteritems(src.tags(iband)) if
                          not key.upper().startswith('STATISTIC')])
                     if len(cleaned_tags) > 0:
                         dest.update_tags(iband, **cleaned_tags)
@@ -1115,8 +1184,10 @@ def multi_block_bands_processing(image_file, pixel_size, out_folder, band_nums=[
 
                     # image statistics have changed so don't copy the tags
                     cleaned_tags = dict(
-                        [(key, val) for key, val in src.tags(iband).iteritems() if
-                         not key.upper().startswith('STATISTIC')])
+                        [(key, val) for key, val in
+                        six.iteritems(src.tags(iband))
+                        if not key.upper().startswith('STATISTIC')]
+                    )
                     if len(cleaned_tags) > 0:
                         dest.update_tags(iband, **cleaned_tags)
 
@@ -1159,8 +1230,11 @@ def multi_block_bands_processing(image_file, pixel_size, out_folder, band_nums=[
                         dest.write(new_band, iband)
 
                         # image statistics have changed so don't copy the tags
-                        cleaned_tags = dict([(key, val) for key, val in src.tags(iband).iteritems()
-                                             if not key.upper().startswith('STATISTIC')])
+                        cleaned_tags = dict([
+                            (key, val) for key, val in
+                            six.iteritems(src.tags(iband))
+                            if not key.upper().startswith('STATISTIC')
+                        ])
 
                         if len(cleaned_tags) > 0:
                             dest.update_tags(iband, **cleaned_tags)
@@ -1183,8 +1257,11 @@ def multi_block_bands_processing(image_file, pixel_size, out_folder, band_nums=[
 
                         dest.write(filled, iband)
                         # image statistics have changed so don't copy the tags
-                        cleaned_tags = dict([(key, val) for key, val in src.tags(iband).iteritems()
-                                             if not key.upper().startswith('STATISTIC')])
+                        cleaned_tags = dict(
+                            [(key, val) for key, val in
+                            six.iteritems(src.tags(iband))
+                            if not key.upper().startswith('STATISTIC')]
+                        )
                         if len(cleaned_tags) > 0:
                             dest.update_tags(iband, **cleaned_tags)
 
@@ -1247,9 +1324,11 @@ def multi_block_bands_processing(image_file, pixel_size, out_folder, band_nums=[
                     dest.write(smooth, 1)
 
                     # image statistics have changed so don't copy the tags
-                    cleaned_tags = dict(
-                        [(key, val) for key, val in src.tags(iband).iteritems() if
-                         not key.upper().startswith('STATISTIC')])
+                    cleaned_tags = dict([
+                        (key, val) for key, val in
+                        six.iteritems(src.tags(iband))
+                        if not key.upper().startswith('STATISTIC')
+                     ])
                     if len(cleaned_tags) > 0:
                         dest.update_tags(1, **cleaned_tags)
 
@@ -1338,7 +1417,7 @@ def calc_indices_for_block(image_file, pixel_size, band_map, out_folder, indices
             raise IOError('{} does not exist-Got {}'.format(ea_arg[0], os.path.dirname(ea_arg[1])))
 
     for ea_arg in [('image_epsg', image_epsg), ('out_epsg', out_epsg)]:
-        if not isinstance(ea_arg[1], (int, long)):
+        if not isinstance(ea_arg[1], six.integer_types):
             raise TypeError('{} must be a integer - Got {}'.format(*ea_arg))
 
     filename, ext = os.path.splitext(os.path.basename(image_file))
@@ -1501,10 +1580,10 @@ def kmeans_clustering(raster_files, output_tif, n_clusters=3, max_iterations=500
         pandas.core.frame.DataFrame: A dataframe containing cluster statistics for each image.
     """
 
-    if not isinstance(n_clusters, (int, long)):
+    if not isinstance(n_clusters, six.integer_types):
         raise TypeError('Size must be an Integer.')
 
-    if not isinstance(max_iterations, (int, long)):
+    if not isinstance(max_iterations, six.integer_types):
         raise TypeError('Size must be an Integer.')
 
     if not isinstance(raster_files, list):
@@ -1657,7 +1736,8 @@ def kmeans_clustering(raster_files, output_tif, n_clusters=3, max_iterations=500
 
     LOGGER.info('{:<30} {:<15} {dur:<15} {}'.format(
         'Saved Stats CSV', '', output_tif.replace('.tif', '_statistics.csv'),
-        dur=timedelta(seconds=time.time() - step_time)))
+        dur=str(timedelta(seconds=time.time() - step_time))
+    ))
 
     # clean up of intermediate files
     if len(temp_file_list) > 0 and not config.get_debug_mode():
@@ -1669,7 +1749,9 @@ def kmeans_clustering(raster_files, output_tif, n_clusters=3, max_iterations=500
 
     LOGGER.info('{:<30}  {:<15} {dur:<15} {}'.format(
         'K-Means Clustering Completed', '', '{} zones for {} rasters'.format(
-            n_clusters, len(raster_files)), dur=timedelta(seconds=time.time() - start_time)))
+            n_clusters, len(raster_files)),
+            dur=str(timedelta(seconds=time.time() - start_time))
+    ))
 
     return results_df
 
@@ -1724,13 +1806,13 @@ def create_points_along_line(lines_geodataframe, lines_crs, distance_between_poi
 
     for argCheck in [('offset_distance', offset_distance),
                      ('distance_between_points', distance_between_points)]:
-        if not isinstance(argCheck[1], (int, long, float)):
+        if not isinstance(argCheck[1], number_types):
             raise TypeError('{} must be a floating number.'.format(argCheck[0]))
 
     if not isinstance(lines_crs, pyprecag_crs.crs):
         raise TypeError('Crs must be an instance of pyprecag.crs.crs')
 
-    if not isinstance(out_epsg, (int, long)):
+    if not isinstance(out_epsg, six.integer_types):
         raise TypeError('out_epsg must be a integer - Got {}'.format(*out_epsg))
 
     if out_points_shapefile is not None and os.path.isabs(out_points_shapefile):
@@ -2055,10 +2137,10 @@ def ttest_analysis(points_geodataframe, points_crs, values_raster, out_folder,
     line_count = len(points_geodataframe['TrialID'].unique())
 
     # Extract raster values for points ------------------------------------------------------------
-
     gdf_points, points_crs = extract_pixel_statistics_for_points(points_geodataframe, points_crs,
                                                                  raster_files, 'extract_pixels.csv',
                                                                  size_list=[1])
+    
     gdf_points.index.name = 'rowID'
     column_names = {}
 
@@ -2084,12 +2166,17 @@ def ttest_analysis(points_geodataframe, points_crs, values_raster, out_folder,
         column_names['Zone'] = 'Zone'
 
     # rename columns from raster name to the dictionary key to be more generic
-    gdf_points.rename(columns=dict([[v, k] for k, v in column_names.iteritems()]), inplace=True)
+    gdf_points.rename(columns=dict([[v, k] for k, v in six.iteritems(column_names)]), inplace=True)
+
+    # group by trialid and find those where all are none ie no overlap with rasters
+    dfnulls = gdf_points.set_index('TrialID')['Value'].isnull().groupby('TrialID').all()
+    
+    # get only the trial id's with data
+    gdf_points = gdf_points[~gdf_points['TrialID'].isin(dfnulls[dfnulls].index.values.tolist())]
+    del dfnulls
 
     # create a unique id for each line/point
-    gdf_points.insert(1, 'TrialPtID',
-                      gdf_points.apply(lambda x: "'{}-{}'".format(x['TrialID'],
-                                                                  x['PointID']), axis=1))
+    gdf_points.insert(1, 'TrialPtID', gdf_points.apply(lambda x: "'{}-{}'".format(x['TrialID'],x['PointID']), axis=1))
 
     # Apply bearing to all transects
     for group_name, gdf_group in gdf_points.groupby(['TrialID', transect_column]):
@@ -2248,9 +2335,9 @@ def ttest_analysis(points_geodataframe, points_crs, values_raster, out_folder,
 
             use_pvalue = 0.05
             df_statstable['sig_color'] = df_statstable['p_value'].apply(
-                lambda x: 'r' if x > use_pvalue else 'k')
+                                                lambda x: 'r' if x > use_pvalue else 'k')
             df_statstable['sig_label'] = df_statstable['p_value'].apply(
-                lambda x: 'not significant' if x > use_pvalue else 'significant')
+                                                lambda x: 'not significant' if x > use_pvalue else 'significant')
 
             # Create the map only once for both lines ----------------------------------------------
             if iscenario == 1:
@@ -2264,10 +2351,10 @@ def ttest_analysis(points_geodataframe, points_crs, values_raster, out_folder,
                 # Aggregate points with the GroupBy - for shapely 1.3.2+ use the point objects
                 try:
                     gdf_lines = gdf_map.groupby(['Strip_Name'])['geometry'].apply(
-                        lambda x: LineString(x.tolist()) if x.size > 1 else x.tolist())
+                                                lambda x: LineString(x.tolist()) if x.size > 1 else x.tolist())
                 except ValueError:
                     gdf_lines = gdf_map.groupby(['Strip_Name'])['geometry'].apply(
-                        lambda x: LineString([(p.x, p.y) for p in x]))
+                                                lambda x: LineString([(p.x, p.y) for p in x]))
 
                 gdf_lines = GeoDataFrame(gdf_lines, geometry='geometry')
 
@@ -2276,10 +2363,8 @@ def ttest_analysis(points_geodataframe, points_crs, values_raster, out_folder,
                 # fix the extent to the vector
                 ax.set_aspect('equal')
                 ax.set_title('Map for Trial {}'.format(line_id), fontsize=8)
-                ax.set_xlim(left=gdf_strip.total_bounds[0] - 50,
-                            right=gdf_strip.total_bounds[2] + 50)
-                ax.set_ylim(bottom=gdf_strip.total_bounds[1] - 50,
-                            top=gdf_strip.total_bounds[3] + 50)
+                ax.set_xlim(left=gdf_strip.total_bounds[0] - 50, right=gdf_strip.total_bounds[2] + 50)
+                ax.set_ylim(bottom=gdf_strip.total_bounds[1] - 50, top=gdf_strip.total_bounds[3] + 50)
 
                 # plot the line
                 gdf_lines.plot(ax=ax, cmap='rainbow', legend=True)
@@ -2290,10 +2375,8 @@ def ttest_analysis(points_geodataframe, points_crs, values_raster, out_folder,
 
                 plt.setp(ax.spines.values(), linewidth=0.5)
 
-                ax.get_xaxis().set_major_formatter(mpl.ticker.FuncFormatter(
-                    lambda x, p: "{:,.0f}".format(x)))
-                ax.get_yaxis().set_major_formatter(mpl.ticker.FuncFormatter(
-                    lambda x, p: "{:,.0f}".format(x)))
+                ax.get_xaxis().set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: "{:,.0f}".format(x)))
+                ax.get_yaxis().set_major_formatter(mpl.ticker.FuncFormatter(lambda x, p: "{:,.0f}".format(x)))
                 # plt.xticks(rotation='vertical')
 
                 anno_kwargs = {'size': 7, 'va': 'center', 'ha': 'center',
diff --git a/pyprecag/pyprecag.py b/pyprecag/pyprecag.py
deleted file mode 100644
index 2678604..0000000
--- a/pyprecag/pyprecag.py
+++ /dev/null
@@ -1,11 +0,0 @@
-# -*- coding: utf-8 -*-
-
-# Ensure backwards compatibility with Python 2
-# from __future__ import (
-#     absolute_import,
-#     division,
-#     print_function,
-#     unicode_literals)
-# from builtins import *
-
-# Need forwards compatible code?
diff --git a/pyprecag/raster_ops.py b/pyprecag/raster_ops.py
index 36a91f0..863f0ca 100644
--- a/pyprecag/raster_ops.py
+++ b/pyprecag/raster_ops.py
@@ -16,6 +16,7 @@
 from rasterio.windows import get_data_window, from_bounds, intersection
 
 from scipy.ndimage import generic_filter
+import six
 from . import crs as pyprecag_crs
 from .bandops import CalculateIndices
 from . import config, TEMPDIR
@@ -47,7 +48,7 @@ def create_raster_transform(bounds, pixel_size, snap_extent_to_pixel=True, buffe
     if not isinstance(bounds, (tuple, list)):
         raise TypeError('Bounds should be a tuple or list for xmin, ymin, xmax, ymax')
 
-    if not isinstance(pixel_size, (int, long, float)):
+    if not isinstance(pixel_size, six.integer_types + (float, )):
         raise TypeError('pixel_size must be numeric number.')
     if not isinstance(snap_extent_to_pixel, (int, bool)):
         raise TypeError('snap_extent_to_pixel must be boolean.')
@@ -99,7 +100,7 @@ def raster_snap_extent(x_min, y_min, x_max, y_max, pixel_size):
 
     for argCheck in [('x_min', x_min), ('y_min', y_min), ('x_max', x_max), ('y_max', y_max),
                      ('pixel_size', pixel_size)]:
-        if not isinstance(argCheck[1], (int, long, float)):
+        if not isinstance(argCheck[1], six.integer_types + (float, )):
             raise TypeError('{} must be a floating number.'.format(argCheck[0]))
 
     b_box = [x_min - (x_min % pixel_size),  # calc new xMin
@@ -139,7 +140,7 @@ def rescale(raster, min_value, max_value, band_num=1, ignore_nodata=True):
                         "rasterio.open(os.path.normpath())")
 
     for argCheck in [('min_value', min_value), ('max_value', max_value)]:
-        if not isinstance(argCheck[1], (int, long, float)):
+        if not isinstance(argCheck[1], six.integer_types + (float, )):
             raise TypeError('{} must be numeric.'.format(argCheck[0]))
 
     band = raster.read(band_num, masked=ignore_nodata)
@@ -292,7 +293,7 @@ def focal_statistics(raster, band_num=1, ignore_nodata=True, size=3, function=np
             warnings.simplefilter("ignore", category=RuntimeWarning)
             filtered = generic_filter(band, function, mode='constant', cval=np.nan, size=size)
 
-        col_name += [function.func_name.replace('nan', ''), '{0}x{0}'.format(size)]
+        col_name += [function.__name__.replace('nan', ''), '{0}x{0}'.format(size)]
     else:
         filtered = band
         col_name += ['pixel']
@@ -447,7 +448,8 @@ def reproject_image(image_file, out_imagefile, out_epsg, band_nums=[],
                               resampling=resampling)
 
                     # image statistics have changed so copy only the other tags
-                    cleaned_tags = dict([(key, val) for key, val in src.tags(i).iteritems()
+                    iter_items = six.iteritems(src.tags(i))
+                    cleaned_tags = dict([(key, val) for key, val in iter_items
                                          if not key.upper().startswith('STATISTIC')])
                     if len(cleaned_tags) > 0:
                         dest.update_tags(i, **cleaned_tags)
diff --git a/pyprecag/table_ops.py b/pyprecag/table_ops.py
index dc2a84b..f53b5dc 100644
--- a/pyprecag/table_ops.py
+++ b/pyprecag/table_ops.py
@@ -1,6 +1,7 @@
 import logging
 
 import pandas as pd
+import six
 from geopandas import GeoDataFrame
 from scipy.stats import stats
 
@@ -70,7 +71,7 @@ def calculate_strip_stats(input_table, treatment_column, control_columns=[], siz
     if not isinstance(control_columns, list):
         raise TypeError("control_columns should be a list.")
 
-    if not isinstance(treatment_column, basestring):
+    if not isinstance(treatment_column, six.string_types):
         raise TypeError("treatment_column should be a string.")
 
     if treatment_column is None or treatment_column == '':
diff --git a/pyprecag/vector_ops.py b/pyprecag/vector_ops.py
index bf16941..4f03a84 100644
--- a/pyprecag/vector_ops.py
+++ b/pyprecag/vector_ops.py
@@ -13,6 +13,7 @@
 import numpy as np
 import pyproj
 import shapely
+import six
 from fiona import collection as fionacoll
 from fiona.crs import to_string
 from osgeo import osr
@@ -49,7 +50,7 @@ def thin_point_by_distance(point_geodataframe, point_crs, thin_distance_metres=1
     """
 
     for argCheck in [('thinDist_m', thin_distance_metres)]:
-        if not isinstance(argCheck[1], (int, long, float)):
+        if not isinstance(argCheck[1], six.integer_types + (float, )):
             raise TypeError('{} must be a floating number.'.format(argCheck[0]))
 
     if not isinstance(point_geodataframe, geopandas.GeoDataFrame):
@@ -153,10 +154,13 @@ def thin_by_distance(curRow, filter_text):
             raise TypeError(
                 "There are no features left after {}. Check the coordinate systems and try again".format(sortBy))
 
-        LOGGER.info(
-            '{:<30} {:>10,}   {dur:<15} {}'.format('Filter by distance - {}'.format(filter_string.replace('point', '')),
-                                                  len(subset), 'del {} pts'.format(stepTotal),
-                                                  dur=datetime.timedelta(seconds=time.time() - filterTime)))
+        LOGGER.info('{:<30} {:>10,}   {dur:<15} {}'.format(
+            'Filter by distance - {}'.format(
+                filter_string.replace('point', '')
+            ),
+            len(subset), 'del {} pts'.format(stepTotal),
+            dur=str(datetime.timedelta(seconds=time.time() - filterTime))
+        ))
 
     # set sort back to original row order
     # point_geodataframe.sort_index(axis=1, ascending=True, inplace=True)
diff --git a/setup.py b/setup.py
index 07ff981..6fa92df 100644
--- a/setup.py
+++ b/setup.py
@@ -61,6 +61,7 @@
         # 'Programming Language :: Python :: 3.2',
         # 'Programming Language :: Python :: 3.3',
         # 'Programming Language :: Python :: 3.4',
+        'Programming Language :: Python :: 3.7',
     ],
 
     # What does your project relate to?
@@ -82,14 +83,15 @@
     # https://packaging.python.org/en/latest/technical.html#install-requires-vs-requirements-files
     install_requires=[
         'future',
+        'scipy',
         'fiona',
-        'gdal', # Required Microsoft Visual C++ Compiler for Python 2.7
-        'rasterio>1,<=1.0.13', # Required Numpy
-        'geopandas>=0.4.0,<0.6.0',
-        'pyproj<2.2.0',
+        'gdal',          # Required Microsoft Visual C++ Compiler for Python 2.7
+        'rasterio',      # Required Numpy
+        'geopandas',
+        'pyproj',
+        'six>1.12.0',
         # These were included to prevent import errors:
         'unidecode',
-        'scipy',
         'chardet',
         'matplotlib',
         ],


Test patch:
diff --git a/pyprecag/tests/make_dummy_data.py b/pyprecag/tests/make_dummy_data.py
index 0f44d14..ffb434c 100644
--- a/pyprecag/tests/make_dummy_data.py
+++ b/pyprecag/tests/make_dummy_data.py
@@ -2,6 +2,27 @@
 import os
 import rasterio
 
+import traceback
+import warnings
+import sys
+
+def warn_with_traceback(message, category, filename, lineno, file=None, line=None):
+    """
+    Add a traceback to warnings to help locate where it being triggered
+    source: https://stackoverflow.com/a/22376126/9567306
+
+    Add the following code
+
+    import warnings
+    warnings.showwarning = make_dummy_data.warn_with_traceback
+
+    Use     warnings.resetwarnings()     to reset
+
+   """
+    log = file if hasattr(file,'write') else sys.stderr
+    traceback.print_stack(file=log)
+    log.write(warnings.formatwarning(message, category, filename, lineno, line))
+
 
 def make_dummy_tif_files(out_dir):
     """ Create fake TIFF to be used as test datasets for other components of pyprecag
diff --git a/pyprecag/tests/test_bandops.py b/pyprecag/tests/test_bandops.py
index dc3239a..a3aa613 100644
--- a/pyprecag/tests/test_bandops.py
+++ b/pyprecag/tests/test_bandops.py
@@ -19,7 +19,7 @@ def setUpClass(cls):
         # 'https://stackoverflow.com/a/34065561'
         super(TestBandOps, cls).setUpClass()
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -44,7 +44,7 @@ def test_band_mapping(self):
         with self.assertRaises(AttributeError) as msg:
             bm = BandMapping(red=4, green=3, blud=2)
         self.assertEqual("BandMapping has no attribute blud. Allowed keys are "
-                "blue, mask, infrared, rededge, green, red", str(msg.exception))
+                "blue, green, infrared, mask, red, rededge", str(msg.exception))
 
         self.assertDictEqual(BandMapping(),
                              {'blue': 0, 'mask': 0, 'infrared': 0, 'rededge': 0, 'green': 0,
diff --git a/pyprecag/tests/test_convert.py b/pyprecag/tests/test_convert.py
index 84b16a9..fb45b83 100644
--- a/pyprecag/tests/test_convert.py
+++ b/pyprecag/tests/test_convert.py
@@ -1,33 +1,35 @@
 import shutil
 import unittest
 
+from pyprecag.tests import make_dummy_data
+
 from pyprecag.convert import *
+from pyprecag import crs as pyprecag_crs
 import time
 import os
 import tempfile
 import logging
 from shapely.geometry import LineString
-
+import geopandas as gpd
 from pyprecag.tests.test_crs import EPSG_28354_WKT
 
 pyFile = os.path.basename(__file__)
-
 TmpDir = tempfile.gettempdir()
 TmpDir = os.path.join(TmpDir, os.path.splitext(pyFile)[0])
+# TmpDir = tempfile.mkdtemp(dir=tempfile.gettempdir(),prefix= os.path.splitext(pyFile)[0] + '_')
 
 this_dir = os.path.abspath(os.path.dirname(__file__))
 
 logging.captureWarnings(True)
 logging.basicConfig(level=logging.INFO, format="%(message)s")
 
-
 class TestConvert(unittest.TestCase):
     @classmethod
     def setUpClass(cls):
         # 'https://stackoverflow.com/a/34065561'
         super(TestConvert, cls).setUpClass()
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -40,6 +42,8 @@ def tearDownClass(cls):
         if not testFailed:
             print ('Tests Passed .. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
+        else:
+            print ('Tests Failed .. {}'.format(TmpDir))
 
     def setUp(self):
         self.startTime = time.time()
@@ -64,7 +68,7 @@ def test_convert_csv_to_points_EastingNorthings(self):
                                                   coord_columns=['Easting', 'Northing'],
                                                   coord_columns_epsg=epsg)
 
-        self.assertIsInstance(gdf_data, GeoDataFrame)
+        self.assertIsInstance(gdf_data, gpd.GeoDataFrame)
         self.assertEqual(len(gdf_data), 13756)
         self.assertEqual(gdf_crs.epsg_number, epsg)
         self.assertTrue(os.path.exists(out_file))
@@ -79,7 +83,7 @@ def test_convert_csv_to_points_WGS84(self):
                                                   out_epsg=out_epsg)
 
         self.assertTrue(os.path.exists(out_file))
-        self.assertIsInstance(gdf_data, GeoDataFrame)
+        self.assertIsInstance(gdf_data,  gpd.GeoDataFrame)
         self.assertEqual(len(gdf_data), 1543)
         self.assertEqual(list(set(gdf_data.geom_type)), ['Point'])
 
@@ -88,7 +92,7 @@ def test_convert_csv_to_points_WGS84(self):
                                        [598868.6709, 6054050.0529, 599242.9049, 6054415.3845], 4)
 
         self.assertEqual(gdf_crs.epsg_number, out_epsg)
-        self.assertEqual(gdf_crs.crs_wkt, EPSG_28354_WKT)
+        self.assertEqual(gdf_crs.crs_wkt[:154], EPSG_28354_WKT[:154])
 
     def test_convert_csv_to_points_WGS84_GuessEPSG(self):
         in_file = os.path.realpath(this_dir + "/data/area2_yield_ISO-8859-1.csv")
@@ -97,11 +101,11 @@ def test_convert_csv_to_points_WGS84_GuessEPSG(self):
         gdf_data, gdf_crs = convert_csv_to_points(in_file, out_file, coord_columns_epsg=4326,
                                                   out_epsg=-1)
         self.assertTrue(os.path.exists(out_file))
-        self.assertIsInstance(gdf_data, GeoDataFrame)
+        self.assertIsInstance(gdf_data,  gpd.GeoDataFrame)
         self.assertEqual(len(gdf_data), 1543)
         self.assertEqual(list(set(gdf_data.geom_type)), ['Point'])
         self.assertEqual(gdf_crs.epsg_number, 28354)
-        self.assertEqual(gdf_crs.crs_wkt, EPSG_28354_WKT)
+        self.assertEqual(gdf_crs.crs_wkt[:154], EPSG_28354_WKT[:154])
 
     def test_numeric_pixelsize_to_string(self):
         self.assertEqual(numeric_pixelsize_to_string(0.42), '42cm')
@@ -125,13 +129,13 @@ def test_cardinal_direction(self):
 
     def test_pts_linebearing(self):
         # line_bearing uses point_to_point_bearing so this doubles as it's tests as well
-        c_line_gdf = {'geometry': [LineString([(740873, 6169764), (741269, 6169764)]),
+        c_line_gdf = gpd.GeoDataFrame({'geometry': [LineString([(740873, 6169764), (741269, 6169764)]),
                                    LineString([(741000, 6169800), (741003, 6170012)]),
                                    LineString([(741401, 6169800), (741057, 6170012)]),
                                    LineString([(740900, 6169912), (740979, 6170094)])],
-                      'LineID': [1, 2, 3, 4]}
+                                    'LineID': [1, 2, 3, 4]}, crs=pyprecag_crs.from_epsg(28354))
 
-        c_line_gdf = GeoDataFrame(c_line_gdf, crs={'init': 'epsg:28354'})
+        c_line_gdf = GeoDataFrame(c_line_gdf, crs=pyprecag_crs.from_epsg(28354))
 
         self.assertEqual(90.0, line_bearing(c_line_gdf.iloc[0]))
         self.assertAlmostEqual(0.8107352192623694, line_bearing(c_line_gdf.iloc[1]), 4)
@@ -139,14 +143,16 @@ def test_pts_linebearing(self):
         self.assertAlmostEqual(23.464022404707464, line_bearing(c_line_gdf.iloc[3]), 4)
 
     def test_drop_z(self):
-        c_line_gdf = {'geometry': [LineString([(740873, 6169764, 5), (741269, 6169764, 5)]),
+        c_line_gdf = gpd.GeoDataFrame( {'geometry': [LineString([(740873, 6169764, 5), (741269, 6169764, 5)]),
                                    LineString([(741000, 6169800), (741003, 6170012)]),
-                                   LineString([(741401.415, 6169800), (741057, 6170012)]),
-                                   LineString([(740900.861, 6169912), (740979, 6170094)])],
-                      'LineID': [1, 2, 3, 4]}
-        c_line_gdf = GeoDataFrame(c_line_gdf, crs={'init': 'epsg:28354'})
-
-
+                                   LineString([(741401.415, 6169800,4), (741057, 6170012,3)]),
+                                   LineString([(740900.861, 6169912,2), (740979, 6170094,5)])],
+                      'LineID': [1, 2, 3, 4]}, crs=pyprecag_crs.from_epsg(28354))
 
+        c_line_gdf['geometry2'] = c_line_gdf['geometry'].apply(lambda x: drop_z(x))
 
+        self.assertEqual(c_line_gdf['geometry'].iloc[1], c_line_gdf['geometry2'].iloc[1])
+        self.assertNotEqual(c_line_gdf['geometry'].iloc[0], c_line_gdf['geometry2'].iloc[0])
 
+        self.assertFalse(c_line_gdf['geometry2'].iloc[0].has_z)
+        self.assertFalse(c_line_gdf['geometry2'].iloc[-1].has_z)
diff --git a/pyprecag/tests/test_crs.py b/pyprecag/tests/test_crs.py
index ee62ee1..6853eeb 100644
--- a/pyprecag/tests/test_crs.py
+++ b/pyprecag/tests/test_crs.py
@@ -4,13 +4,17 @@
 import time
 import shutil
 import tempfile
+
+import pyproj
 import rasterio
 import logging
 
+from pkg_resources import parse_version
+
 from pyprecag.tests import make_dummy_data
-from pyprecag.crs import crs, getProjectedCRSForXY, getCRSfromRasterFile, getUTMfromWGS84, \
-    distance_metres_to_dd
-from fiona.crs import from_epsg
+from pyprecag.crs import crs, getProjectedCRSForXY, getCRSfromRasterFile, getUTMfromWGS84, distance_metres_to_dd
+
+from pyprecag.crs import from_epsg
 
 this_dir = os.path.abspath(os.path.dirname(__file__))
 
@@ -73,7 +77,7 @@ def setUpClass(cls):
         super(TestCrsClass, cls).setUpClass()
 
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -99,11 +103,10 @@ def test_getFromEPSG(self):
         test.getFromEPSG(28354)
         self.assertEqual(test.epsg_number, 28354)
         self.assertEqual(test.epsg, from_epsg(test.epsg_number))
-        self.assertEqual(EPSG_28354_WKT, test.crs_wkt)
+        self.assertEqual(EPSG_28354_WKT[:154], test.crs_wkt[:154])
 
-        self.assertEqual('+proj=utm +zone=54 +south +ellps=GRS80 '
-                         '+towgs84=0,0,0,0,0,0,0 +units=m +no_defs ',
-                         test.proj4)
+        self.assertEqual('+proj=utm +zone=54 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs',
+                         test.proj4.strip())
 
     def test_getFromWKT_GDA1(self):
         test = crs()
@@ -111,75 +114,72 @@ def test_getFromWKT_GDA1(self):
 
         self.assertEqual(test.epsg_number, 28354)
         self.assertEqual(test.epsg, from_epsg(test.epsg_number))
-        self.assertEqual(test.proj4, '+proj=utm +zone=54 +south +ellps=GRS80 '
-                                     '+towgs84=0,0,0,0,0,0,0 +units=m +no_defs ')
+        self.assertEqual('+proj=utm +zone=54 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs',
+                         test.proj4.strip())
+
         self.assertEqual(EPSG_28354_WKT, test.crs_wkt)
 
         test.getFromWKT(ESRI_54_WKT_2)
         self.assertEqual(28354, test.epsg_number)
         self.assertEqual(from_epsg(test.epsg_number), test.epsg)
-        self.assertEqual('+proj=utm +zone=54 +south +ellps=GRS80 +units=m +no_defs ', test.proj4)
+        self.assertEqual('+proj=utm +zone=54 +south +ellps=GRS80 +units=m +no_defs', test.proj4.strip())
         self.assertEqual(ESRI_54_WKT_2, test.crs_wkt)
 
     def test_getUTMfromWGS84(self):
-        result = getUTMfromWGS84(138.679870, -34.037740)
-        self.assertEqual(54, result[0])
-        self.assertEqual(ESRI_54_SUTM_WKT, result[1].ExportToWkt())
 
-    def test_distance_metres_to_dd(self):
-        dist = distance_metres_to_dd(138.822027994089, -34.4842175261199, 500)
-        self.assertEqual(round(dist, 8), 0.00544233)
+        utm_zone, _, _, out_epsg = getUTMfromWGS84(138.679870, -34.037740)
 
-    def testGetFromWKTUsingLookup_GDA(self):
+        self.assertEqual(54, utm_zone)
+        self.assertEqual(32754, out_epsg)
 
-        with warnings.catch_warnings(record=True) as w:
-            test = crs()
-            test.getFromWKT(EPSG_28354_WKT)
+        utm_zone, _, _, out_epsg = getUTMfromWGS84(138.679870, 34.037740)
 
-            if len(w) == 1:
-                print ('Timeout Occurred')
-                self.assertIn('could not be reached. Timeout after 10 seconds', str(w[-1].message))
-            else:
-                self.assertEqual(EPSG_28354_WKT, test.crs_wkt)
+        self.assertEqual(54, utm_zone)
+        self.assertEqual(32654, out_epsg)
 
-                self.assertEqual(test.epsg, from_epsg(test.epsg_number))
-                # self.assertEqual(test.proj4, '+proj=utm +zone=54 +south +ellps=GRS80'
-                #                              ' +units=m +no_defs ')
-                self.assertEqual(test.epsg_predicted, False)
+        # Don't test wkt's anymore, they keep changing minutely
+        #self.assertEqual(ESRI_54_SUTM_WKT, result[1].ExportToWkt())
 
-            test = crs()
-            test.getFromWKT(ESRI_54_WKT_1)
+    def test_distance_metres_to_dd(self):
+        dist = distance_metres_to_dd(138.822027994089, -34.4842175261199, 500)
+        self.assertAlmostEqual(0.00544233, dist, 8)
 
-            if len(w) == 1:
-                print ('Timeout Occurred')
-                self.assertIn('could not be reached. Timeout after 10 seconds', str(w[-1].message))
-            else:
-                self.assertEqual(ESRI_54_WKT_1, test.crs_wkt)
+    def testGetFromWKTUsingLookup_GDA(self):
+        test = crs()
+        test.getFromWKT(EPSG_28354_WKT)
 
-                self.assertEqual(test.epsg, from_epsg(test.epsg_number))
-                # self.assertEqual(test.proj4, '+proj=utm +zone=54 +south +ellps=GRS80'
-                #                              ' +units=m +no_defs ')
-                self.assertEqual(test.epsg_predicted, True)
+        self.assertEqual(EPSG_28354_WKT, test.crs_wkt)
+        self.assertEqual(test.epsg, from_epsg(test.epsg_number))
+        # self.assertEqual(test.proj4, '+proj=utm +zone=54 +south +ellps=GRS80 +units=m +no_defs ')
+        self.assertEqual(test.epsg_predicted, False)
 
+        test = crs()
+        test.getFromWKT(ESRI_54_WKT_1)
 
-            test = crs()
-            test.getFromWKT(ESRI_54_WKT_2)
+        self.assertEqual(ESRI_54_WKT_1, test.crs_wkt)
+        self.assertEqual(test.epsg, from_epsg(test.epsg_number))
+        # self.assertEqual(test.proj4, '+proj=utm +zone=54 +south +ellps=GRS80 +units=m +no_defs ')
+        if parse_version(pyproj.__version__) >= parse_version('2.5.0'):
+            self.assertEqual(test.epsg_predicted, False)
+        else:
+            self.assertEqual(test.epsg_predicted, True)
 
-            if len(w) == 1:
-                print ('Timeout Occurred')
-                self.assertIn('could not be reached. Timeout after 10 seconds', str(w[-1].message))
-            else:
-                self.assertEqual(ESRI_54_WKT_2, test.crs_wkt)
+        test = crs()
+        test.getFromWKT(ESRI_54_WKT_2)
 
-                self.assertEqual(test.epsg, from_epsg(test.epsg_number))
-                # self.assertEqual(test.proj4, '+proj=utm +zone=54 +south +ellps=GRS80'
-                #                              ' +units=m +no_defs ')
-                self.assertEqual(test.epsg_predicted, True)
+        self.assertEqual(ESRI_54_WKT_2, test.crs_wkt)
+        self.assertEqual(test.epsg, from_epsg(test.epsg_number))
+        # self.assertEqual(test.proj4, '+proj=utm +zone=54 +south +ellps=GRS80 +units=m +no_defs ')
+        if parse_version(pyproj.__version__) >= parse_version('2.5.0'):
+            self.assertEqual(test.epsg_predicted, False)
+        else:
+            self.assertEqual(test.epsg_predicted, True)
 
     def test_getProjectedCRSForXY(self):
+
         # wgs84 edge of z54-55
         result = getProjectedCRSForXY(143.95231, -37.79412, 4326)
-        self.assertEqual(EPSG_28354_WKT, result.srs.ExportToWkt())
+        self.assertEqual(EPSG_28354_WKT[:154], result.srs.ExportToWkt()[:154])
 
         # Same point as above but in AGD66
         result = getProjectedCRSForXY(143.95099, -37.79561, 4202)
@@ -188,14 +188,11 @@ def test_getProjectedCRSForXY(self):
         # a New Zealand point
         result = getProjectedCRSForXY(169.796934, -44.380541)
         self.assertEqual(result.epsg_number, 2193)
-        self.assertEqual(NZ_WKT, result.crs_wkt)
+        self.assertEqual(NZ_WKT[:183], result.crs_wkt[:183])
 
         # A northern hemisphere point
         result = getProjectedCRSForXY(143.95099, 37.79561, 4326)
         self.assertEqual(result.epsg_number, 32654)
-        self.assertIn('PROJCS["UTM Zone 54, Northern Hemisphere",GEOGCS["WGS 84",'
-                      'DATUM["WGS_1984",SPHEROID["WGS 84"',
-                      result.crs_wkt)
 
         # a southern hemisphere point outside australia
         result = getProjectedCRSForXY(165.95099, -37.79561, 4326)
diff --git a/pyprecag/tests/test_describe.py b/pyprecag/tests/test_describe.py
index 364411d..0dcee97 100644
--- a/pyprecag/tests/test_describe.py
+++ b/pyprecag/tests/test_describe.py
@@ -5,10 +5,8 @@
 from pkg_resources import parse_version, get_distribution
 
 from pyprecag.describe import predictCoordinateColumnNames, CsvDescribe, VectorDescribe
-
 this_dir = os.path.abspath(os.path.dirname(__file__))
 
-
 class test_GeneralDescribe(TestCase):
     def setUp(self):
         self.startTime = time.time()
@@ -36,6 +34,8 @@ def test_from_ISO_8859_1_csv(self):
         descCSV = CsvDescribe(file)
         self.assertEqual(predictCoordinateColumnNames(descCSV.get_column_names()), ['Longitude', 'Latitude'])
 
+        df = descCSV.open_pandas_dataframe()
+        self.assertListEqual(df.columns.to_list(), descCSV.get_column_names())
 
 class test_VectorDescribe_QGIS(TestCase):
     def setUp(self):
@@ -78,9 +78,35 @@ def test_wgs84_mixedPartLine_MZ_esriprj(self):
         self.assertTrue(vDesc.is_mz_aware)
         self.assertEqual(vDesc.geometry_type, 'MultiLineString')
         self.assertEqual(vDesc.feature_count,8)
+        from collections import OrderedDict
+
+        compare_to = OrderedDict([(u'Id', {'shapefile': u'Id', 'alias': 'Id', 'type': 'int', 'dtype': 'int64'}),
+                      (u'float_3dLe', {'shapefile': u'float_3dLe', 'alias': 'float_3dLe', 'type': 'float', 'dtype': 'float64'}),
+                      (u'double_Len', {'shapefile': u'double_Len', 'alias': 'double_Len', 'type': 'float', 'dtype': 'float64'}),
+                     (u'short_id', {'shapefile': u'short_id', 'alias': 'short_id', 'type': 'int', 'dtype': 'int64'}),
+                     (u'Long_vert', {'shapefile': u'Long_vert', 'alias': 'Long_vert', 'type': 'int', 'dtype': 'int64'}),
+                     (u'Date_Creat',{'shapefile': u'Date_Creat', 'alias': 'Date_Creat', 'type': 'str', 'dtype': 'object'}),
+                     (u'part_type', {'shapefile': u'part_type', 'alias': 'part_type', 'type': 'str', 'dtype': 'object'}),
+                     ('geometry',{'shapefile': 'geometry', 'alias': 'geometry', 'type': 'geometry', 'dtype': 'object'})])
+
+        self.assertListEqual(list(vDesc.column_properties.keys()), list(compare_to.keys()))
+
+        import geopandas
+        for key,val in vDesc.column_properties.items():
+            if hasattr(geopandas, 'array') and key == 'geometry':
+                compare_to['geometry']['dtype'] = 'geometry'
+
+            self.assertDictEqual(vDesc.column_properties[key], val)
+
+    def test_wgs84_mixedPartLine_MZ_esriprj_GH(self):
+        vDesc = VectorDescribe(os.path.realpath(this_dir + "/data/LineMZ_wgs84_MixedPartFieldsTypes_esri.shp"))
+        self.assertEqual(vDesc.crs.epsg_number, 4326)
+        self.assertTrue(vDesc.is_mz_aware)
+        self.assertEqual(vDesc.geometry_type, 'MultiLineString')
+        self.assertEqual(vDesc.feature_count,8)
 
         from collections import OrderedDict
-        self.maxDiff
+        self.maxDiff =None
         if parse_version(get_distribution('geopandas').version) > parse_version('0.6.0'):
             self.assertDictEqual(vDesc.column_properties, OrderedDict(
                 [(u'Id', {'shapefile': u'Id', 'alias': 'Id', 'type': 'int', 'dtype': 'int64'}),
@@ -103,7 +129,7 @@ def test_wgs84_mixedPartLine_MZ_esriprj(self):
                                                                   (u'Date_Creat', {'shapefile': u'Date_Creat', 'alias': 'Date_Creat', 'type': 'str', 'dtype': 'object'}),
                                                                   (u'part_type', {'shapefile': u'part_type', 'alias': 'part_type', 'type': 'str', 'dtype': 'object'}),
                                                                   ('geometry', {'shapefile': 'geometry', 'alias': 'geometry', 'type': 'geometry', 'dtype': 'object'})]))
-
+                                                                  
     def test_wgs84_mixedPartPoly_MZ_esriprj(self):
         vDesc = VectorDescribe(os.path.realpath(this_dir + "/data/PolyMZ_wgs84_MixedPartFieldsTypes.shp"))
         self.assertEqual(vDesc.crs.epsg_number, 4326)
@@ -122,12 +148,28 @@ def tearDown(self):
     def test_csvfile_UTF8(self):
         csvDesc = CsvDescribe(os.path.realpath(this_dir + "/data/area2_yield_ISO-8859-1.csv"))
 
-        self.assertEqual(csvDesc.file_encoding, 'ISO-8859-1')
+
+        # chardet seems to detect this file as ISO-8869-9, which is described
+        # as "Largely the same as ISO/IEC 8859-1, replacing the rarely used
+        # Icelandic letters with Turkish ones."
+
+        #self.assertEqual(csvDesc.file_encoding, 'ISO-8859-1')
+        self.assertEqual(csvDesc.file_encoding, 'ISO-8859-9')
         self.assertEqual(csvDesc.row_count, 1543)
         self.assertEqual(csvDesc.column_count, 18)
         self.assertEqual(predictCoordinateColumnNames(csvDesc.get_column_names()), ['Longitude', 'Latitude'])
         self.assertTrue(csvDesc.has_column_header)
 
+        self.assertEqual(csvDesc.get_column_names()[-1], csvDesc.get_alias_column_names()[-1])
+        self.assertNotEqual(csvDesc.get_column_names()[-2], csvDesc.get_alias_column_names()[-2])
+
+        self.assertEqual(u'Crop Flw(V)(m\xb3/s)', csvDesc.get_column_names()[-2])
+        self.assertEqual('CropFlw(V)(m3/s)', csvDesc.get_alias_column_names()[-2])
+
+        #check to see if unicode characters exist True if all ascii, false if not
+        self.assertTrue(all(ord(char) < 128 for char in csvDesc.get_alias_column_names()[-2]))
+        self.assertFalse(all(ord(char) < 128 for char in csvDesc.get_column_names()[-2]))
+
     def test_csvfile_ascii(self):
         csvDesc = CsvDescribe(os.path.realpath(this_dir + "/data/area1_yield_ascii_wgs84.csv"))
 
diff --git a/pyprecag/tests/test_end2end.py b/pyprecag/tests/test_end2end.py
index 5398c33..a0b881f 100644
--- a/pyprecag/tests/test_end2end.py
+++ b/pyprecag/tests/test_end2end.py
@@ -9,7 +9,7 @@
 from pyprecag.describe import CsvDescribe, predictCoordinateColumnNames
 from pyprecag.processing import *
 from pyprecag.raster_ops import rescale, normalise
-from pyprecag.kriging_ops import prepare_for_vesper_krige, vesper_text_to_raster, run_vesper
+from pyprecag.kriging_ops import prepare_for_vesper_krige, vesper_text_to_raster, run_vesper, VesperControl
 
 pyFile = os.path.basename(__file__)
 
@@ -40,7 +40,7 @@ def setUpClass(cls):
         # 'https://stackoverflow.com/a/34065561'
         super(TestEnd2End, cls).setUpClass()
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -115,16 +115,17 @@ def test03_vectorDescribe(self):
     def test04_blockGrid(self):
 
         global fileBlockTif, file_block_txt
-        fileBlockTif = os.path.join(TmpDir, os.path.splitext(
-            os.path.basename(file_csv))[0] + '_block.tif')
-        file_block_txt = os.path.join(TmpDir, os.path.splitext(
-            os.path.basename(file_csv))[0] + '_block_v.txt')
+        fileBlockTif = os.path.join(TmpDir, os.path.splitext(os.path.basename(file_csv))[0] + '_block.tif')
+        file_block_txt = os.path.join(TmpDir, os.path.splitext(os.path.basename(file_csv))[0] + '_block_v.txt')
+
+        vect_desc = VectorDescribe(fileBox)
 
         if not os.path.exists(fileBlockTif):
             block_grid(in_shapefilename=fileBox,
                        pixel_size=2.5,
                        out_rasterfilename=fileBlockTif,
                        out_vesperfilename=file_block_txt,
+                       out_epsg=vect_desc.crs.epsg_number,
                        snap=True,
                        overwrite=True)
 
@@ -158,6 +159,7 @@ def test05_cleanTrimPoints(self):
         gdf_points, gdfpts_crs = convert.convert_csv_to_points(file_csv,
                                                                coord_columns_epsg=4326,
                                                                out_epsg=epsg)
+
         gdf_out, crs_out = clean_trim_points(gdf_points, gdfpts_crs, data_col, fileTrimmed,
                                              out_keep_shapefile=file_shp,
                                              out_removed_shapefile=file_removed,
@@ -167,7 +169,7 @@ def test05_cleanTrimPoints(self):
         self.assertTrue(os.path.exists(fileTrimmed))
         self.assertTrue(os.path.exists(file_shp))
         self.assertTrue(os.path.exists(file_removed))
-        self.assertEqual(gdf_out.crs, {'init': 'epsg:{}'.format(epsg), 'no_defs': True})
+        self.assertEqual(gdf_out.crs, crs.from_epsg(epsg))
         self.assertEqual(len(gdf_out), 648)
         self.assertIn('nrm_' + data_col, gdf_out.columns)
         self.assertIn('Easting', gdf_out.columns)
@@ -186,11 +188,15 @@ def test06_prepareForVesperKrig(self):
         sub_file = os.path.splitext(os.path.basename(file_csv))[0]
         file_control = sub_file + '_control_' + data_col + '.txt'
 
+        vc = VesperControl()
+        vc.update(xside=30, yside=30)
+
         if not os.path.exists(file_control):
             bat_file, file_control = prepare_for_vesper_krige(df_csv, data_col, file_block_txt,
-                                                              TmpDir, block_size=30,
+                                                              TmpDir,
                                                               control_textfile=file_control,
-                                                              coord_columns=[], epsg=epsg)
+                                                              coord_columns=[], epsg=epsg,
+                                                              control_options=vc)
 
             self.assertTrue(os.path.exists(bat_file))
             self.assertTrue(os.path.exists(file_control))
@@ -309,14 +315,14 @@ def test11_rescaleNormaliseRaster(self):
         with rasterio.open(os.path.normpath(out_normalised), 'w', **out_meta) as out:
             out.write_band(1, rescaled2)
 
-        self.assertAlmostEqual(float(np.nanmax(norm)), 2.0000722408294678, 4)
-        self.assertAlmostEqual(float(np.nanmin(norm)), -2.266947031021118, 4)
+        self.assertAlmostEqual(2.0000722408294678, float(np.nanmax(norm)), 4)
+        self.assertAlmostEqual(-2.266947031021118, float(np.nanmin(norm)), 4)
 
-        self.assertEqual(np.nanmin(rescaled), 0)
-        self.assertEqual(np.nanmax(rescaled), 255)
+        self.assertEqual(0, np.nanmin(rescaled) )
+        self.assertEqual(255, np.nanmax(rescaled))
 
-        self.assertEqual(np.nanmin(rescaled2), 0)
-        self.assertEqual(np.nanmax(rescaled2), 5)
+        self.assertEqual(0, np.nanmin(rescaled2))
+        self.assertEqual(5, np.nanmax(rescaled2))
 
     def test12_kmeansCluster(self):
         out_img = os.path.join(TmpDir, 'kmeans-cluster_3cluster_3rasters.tif')
@@ -340,10 +346,16 @@ def test12_kmeansCluster(self):
 
         with rasterio.open(out_img) as src:
             self.assertEqual(1, src.count)
+            if hasattr(src.crs, 'to_proj4'):
+                self.assertEqual(src.crs.to_proj4(), '+init=epsg:28354')
+            else:
+                self.assertEqual(src.crs.to_string(), '+init=epsg:28354')
             self.assertIn(src.crs.to_string(), ['EPSG:28354', '+init=epsg:28354'])
             self.assertEqual(0, src.nodata)
             band1 = src.read(1, masked=True)
-            self.assertItemsEqual(np.array([0, 1, 2, 3]), np.unique(band1.data))
+
+            six.assertCountEqual(self, np.array([0, 1, 2, 3]), np.unique(band1.data))
+
 
     def test99_gridExtract(self):
         out_fold = os.path.join(TmpDir, 'gridextract')
diff --git a/pyprecag/tests/test_kriging_ops.py b/pyprecag/tests/test_kriging_ops.py
index 384180c..fa38ee4 100644
--- a/pyprecag/tests/test_kriging_ops.py
+++ b/pyprecag/tests/test_kriging_ops.py
@@ -9,6 +9,8 @@
 import rasterio
 import time
 
+from pyprecag.tests import make_dummy_data
+
 from pyprecag import convert, processing
 from pyprecag.describe import predictCoordinateColumnNames, CsvDescribe
 from pyprecag.processing import clean_trim_points
@@ -25,6 +27,7 @@
 logging.captureWarnings(True)
 logging.basicConfig(level=logging.INFO, format="%(message)s")
 
+
 class TestVesperControl(unittest.TestCase):
     def test_updated_keys(self):
         vc = VesperControl(title='This is a test tile')
@@ -110,7 +113,7 @@ def setUpClass(cls):
         super(TestKrigingOps, cls).setUpClass()
 
         if os.path.exists(TEMPDIR):
-            print 'Folder Exists.. Deleting {}'.format(TEMPDIR)
+            print('Folder Exists.. Deleting {}'.format(TEMPDIR))
             shutil.rmtree(TEMPDIR)
 
         os.mkdir(TEMPDIR)
@@ -137,33 +140,6 @@ def run(self, result=None):
         if len(result.failures) > 0 or len(result.errors) > 0:
             test_failed = True
 
-    def test1a_CreateControlHighDensity_block_size(self):
-        file_csv = os.path.realpath(this_dir + "/area2_high_trimmed.csv")
-        grid_file = os.path.realpath(this_dir + "/rasters/area2_5m_blockgrid_v.txt")
-        data_col = r'Yld Mass(Dry)(tonne/ha)'
-
-        csv_desc = CsvDescribe(file_csv)
-        df_csv = csv_desc.open_pandas_dataframe()
-
-        # check using block_size argument - backwards compatible
-        file_bat, file_ctrl = prepare_for_vesper_krige(df_csv, data_col, grid_file, TEMPDIR,
-                                                       control_textfile='test_high_5m_control.txt',
-                                                       block_size=30, coord_columns=[],
-                                                       epsg=28354)
-        if os.path.exists(kriging_ops.vesper_exe):
-            self.assertTrue(os.path.exists(os.path.join(TEMPDIR, 'Vesper/Do_Vesper.bat')))
-        else:
-            self.assertEqual('', file_bat)
-
-        self.assertTrue(os.path.exists(os.path.join(TEMPDIR, 'Vesper', 'test_high_5m_control.txt')))
-        self.assertTrue(os.path.exists(os.path.join(TEMPDIR, 'Vesper', 'test_high_5m_vesperdata.csv')))
-
-        src_df = pd.read_csv(os.path.realpath(this_dir +'/VESPER/high_5m_vesperdata.csv'))
-        test_df = pd.read_csv(os.path.join(TEMPDIR, 'Vesper', 'test_high_5m_vesperdata.csv'))
-
-        pd.testing.assert_frame_equal( src_df,test_df)
-
-
     def test1_CreateControlHighDensity_VesperControlClass(self):
         # check using VesperControl class
         file_csv = os.path.realpath(this_dir + "/area2_high_trimmed.csv")
@@ -211,7 +187,7 @@ def test2_prepareForVesperKrig_LowDensity(self):
                           'iwei': 'no_pairs/variance', 'CO': 92.71, 'C1': 277.9, 'A1': 116.0 })
 
         file_bat, file_ctrl = prepare_for_vesper_krige(csv_desc.open_pandas_dataframe(),
-                                                       data_col, grid_file , TEMPDIR,
+                                                       data_col, grid_file, TEMPDIR,
                                                        control_textfile='test_low_control.txt',
                                                        control_options=ctrl_para,
                                                        coord_columns=[], epsg=28354)
diff --git a/pyprecag/tests/test_processing.py b/pyprecag/tests/test_processing.py
index 6003b36..1ddc7d6 100644
--- a/pyprecag/tests/test_processing.py
+++ b/pyprecag/tests/test_processing.py
@@ -2,8 +2,6 @@
 import tempfile
 import unittest
 
-from geopandas import GeoSeries, GeoDataFrame
-
 from pyprecag.bandops import BandMapping, CalculateIndices
 from pyprecag.crs import crs
 from pyprecag.tests import make_dummy_data
@@ -16,7 +14,7 @@
 
 this_dir = os.path.abspath(os.path.dirname(__file__))
 
-logging.captureWarnings(False)
+logging.captureWarnings(True)
 logging.basicConfig(level=logging.INFO, format="%(message)s")
 
 
@@ -27,7 +25,7 @@ def setUpClass(cls):
         super(TestProcessing, cls).setUpClass()
 
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -58,16 +56,18 @@ def test_BlockGrid(self):
         poly = os.path.realpath(this_dir + "/data/area2_onebox_94mga54.shp")
 
         file_sub_name = os.path.join(TmpDir, os.path.splitext(os.path.basename(poly))[0])
-
-        block_grid(in_shapefilename=poly,
+        vect_desc = VectorDescribe(poly)
+        output_files = block_grid(in_shapefilename=poly,
                    pixel_size=5,
                    out_rasterfilename=file_sub_name + '_block.tif',
                    out_vesperfilename=file_sub_name + '_block_v.txt',
+                   out_epsg=vect_desc.crs.epsg_number,
                    snap=True,
                    overwrite=True)
 
         self.assertTrue(os.path.exists(file_sub_name + '_block.tif'))
         self.assertTrue(os.path.exists(file_sub_name + '_block_v.txt', ))
+        self.assertTrue(len(output_files),1)
 
         with rasterio.open(os.path.normpath(file_sub_name + '_block.tif')) as dataset:
             self.assertEqual(dataset.count, 1)
@@ -76,6 +76,31 @@ def test_BlockGrid(self):
             self.assertEqual(dataset.nodatavals, (-9999.0,))
             self.assertEqual(dataset.dtypes, ('int16',))
 
+    def test_BlockGrid_GrpBy(self):
+        poly = os.path.realpath(this_dir + "/data/PolyMZ_wgs84_MixedPartFieldsTypes.shp")
+
+        file_sub_name = os.path.join(TmpDir, os.path.splitext(os.path.basename(poly))[0])
+
+        output_files = block_grid(in_shapefilename=poly,
+                   pixel_size=5,
+                   out_rasterfilename=file_sub_name + '_block.tif',
+                   out_vesperfilename=file_sub_name + '_block_v.txt',
+                   out_epsg=28354,
+                   groupby='part_type',
+                   snap=True,
+                   overwrite=True)
+
+        self.assertTrue(len(output_files), 2)
+        self.assertTrue(os.path.exists(file_sub_name + '_block_MultiPart.tif'))
+        self.assertTrue(os.path.exists(file_sub_name + '_block_SinglePart_v.txt', ))
+
+        with rasterio.open(os.path.normpath(file_sub_name + '_block_MultiPart.tif')) as dataset:
+            self.assertEqual(dataset.count, 1)
+            self.assertEqual(dataset.width, 59)
+            self.assertEqual(dataset.height, 28)
+            self.assertEqual(dataset.nodatavals, (-9999.0,))
+            self.assertEqual(dataset.dtypes, ('int16',))
+
     def test_cleanTrimPoints(self):
         in_csv = os.path.join(this_dir + "/data/area2_yield_ISO-8859-1.csv")
         in_poly = os.path.join(this_dir + "/data/area2_onebox_94mga54.shp")
@@ -93,7 +118,7 @@ def test_cleanTrimPoints(self):
         self.assertIsInstance(out_gdf, GeoDataFrame)
         self.assertTrue(os.path.exists(out_csv))
         self.assertTrue(gdf_pts_crs, out_crs)
-        self.assertEqual(out_gdf.crs, {'init': 'epsg:28354', 'no_defs': True})
+        self.assertEqual(out_gdf.crs, from_epsg(28354)) #  {'init': 'epsg:28354', 'no_defs': True})
         self.assertEqual(len(out_gdf), 554)
         self.assertIn('EN_EPSG', out_gdf.columns)
 
@@ -158,11 +183,15 @@ def test_kmeansCluster(self):
 
         with rasterio.open(out_img) as src:
             self.assertEqual(1, src.count)
-            self.assertIn(src.crs.to_string(), ['EPSG:28354','+init=epsg:28354'])
-
+            if hasattr(src.crs,'to_proj4'):
+                self.assertEqual(src.crs.to_proj4(), '+init=epsg:28354')
+            else:
+                self.assertEqual(src.crs.to_string(), '+init=epsg:28354')
             self.assertEqual(0, src.nodata)
             band1 = src.read(1, masked=True)
-            self.assertItemsEqual(np.array([0, 1, 2, 3]), np.unique(band1.data))
+            six.assertCountEqual(
+                    self, np.array([0, 1, 2, 3]), np.unique(band1.data)
+            )
 
     def test_PersistorAllYears(self):
         raster_files = glob.glob(os.path.realpath(this_dir + '/data/rasters/Year*.tif'))
@@ -189,7 +218,7 @@ def test_PersistorAllYears(self):
             np.testing.assert_array_equal(src.read(), test.read())
 
             band1 = test.read(1, masked=True)
-            self.assertItemsEqual([-9999, 0, 1, 2, 3], np.unique(band1.data).tolist())
+            six.assertCountEqual(self, [-9999, 0, 1, 2, 3], np.unique(band1.data).tolist())
 
     def test_PersistorTargetProb(self):
         raster_files = glob.glob(os.path.realpath(this_dir + '/data/rasters/Year*.tif'))
@@ -209,7 +238,7 @@ def test_PersistorTargetProb(self):
             np.testing.assert_array_equal(src.read(), test.read())
 
             band1 = test.read(1, masked=True)
-            self.assertItemsEqual([-9999, -1, 0, 1], np.unique(band1.data).tolist())
+            six.assertCountEqual(self, [-9999, -1, 0, 1], np.unique(band1.data).tolist())
 
 
 class TestStripTrials(unittest.TestCase):
@@ -219,7 +248,7 @@ def setUpClass(cls):
         super(TestStripTrials, cls).setUpClass()
 
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -262,7 +291,7 @@ def test_CreateStripTreatmentPoints(self):
                                       (741300, 6169800, 7), (741250, 6169950, 7),
                                       (741150, 6169950, 7), (741100, 6170000, 7)]),
                           LineString([(740800, 6169912, 8), (740900, 6170094, 8)])]
-                , 'TrialID': [1, 2, 3, 4]}, crs={'init': 'epsg:28354'})
+                , 'TrialID': [1, 2, 3, 4]}, crs=pyprecag_crs.from_epsg(28354))
 
         gdf_lines_crs = crs()
         gdf_lines_crs.getFromEPSG(28354)
@@ -368,7 +397,7 @@ def test_TTestAnalysis(self):
                 [127, 3, "SE Strip", 13, 247.6, Point(300806, 6181771)]]
 
         gdf_points = GeoDataFrame(data, columns=columns, geometry='geometry',
-                                  crs={'init': 'epsg:28354'})
+                                  crs=from_epsg(28354))
         crs_points = crs()
         crs_points.getFromEPSG(28354)
 
@@ -398,7 +427,7 @@ def setUpClass(cls):
         super(TestExtractRasterStatisticsForPoints, cls).setUpClass()
 
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -479,7 +508,7 @@ def setUpClass(cls):
         # 'https://stackoverflow.com/a/34065561'
         super(TestCalculateImageIndices, cls).setUpClass()
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.makedirs(TmpDir)
@@ -638,7 +667,7 @@ def setUpClass(cls):
         super(TestResampleToBlock, cls).setUpClass()
 
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.makedirs(TmpDir)
diff --git a/pyprecag/tests/test_rasterops.py b/pyprecag/tests/test_rasterops.py
index 3a97a71..cb0e094 100644
--- a/pyprecag/tests/test_rasterops.py
+++ b/pyprecag/tests/test_rasterops.py
@@ -4,6 +4,7 @@
 import time
 import unittest
 import rasterio
+import six
 
 from pyprecag.tests import make_dummy_data
 from pyprecag import raster_ops
@@ -21,7 +22,7 @@ def setUpClass(cls):
         # 'https://stackoverflow.com/a/34065561'
         super(test_rasterOps, cls).setUpClass()
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
@@ -47,8 +48,14 @@ def test_rasterSnapExtent(self):
         xmin, ymin = 350127.023547, 6059756.84652457
         xmax, ymax = xmin + 3.142 * 200, ymin + 3.142 * 550
 
-        self.assertItemsEqual([350125.0, 6059755.0, 350760.0, 6061485.0], raster_ops.raster_snap_extent(xmin, ymin, xmax, ymax, 5))
-        self.assertItemsEqual([350127.0, 6059756.5, 350755.5, 6061485.0], raster_ops.raster_snap_extent(xmin, ymin, xmax, ymax, 0.5))
+        six.assertCountEqual(
+            self, [350125.0, 6059755.0, 350760.0, 6061485.0],
+            raster_ops.raster_snap_extent(xmin, ymin, xmax, ymax, 5)
+        )
+        six.assertCountEqual(
+            self, [350127.0, 6059756.5, 350755.5, 6061485.0],
+            raster_ops.raster_snap_extent(xmin, ymin, xmax, ymax, 0.5)
+        )
 
     def test_rescaleSingleBand(self):
         with rasterio.open(os.path.normpath(self.singletif)) as src:
diff --git a/pyprecag/tests/test_vector_ops.py b/pyprecag/tests/test_vector_ops.py
index 74c2fc2..b556181 100644
--- a/pyprecag/tests/test_vector_ops.py
+++ b/pyprecag/tests/test_vector_ops.py
@@ -22,7 +22,7 @@ def setUpClass(cls):
         super(test_vectorOps, cls).setUpClass()
 
         if os.path.exists(TmpDir):
-            print 'Folder Exists.. Deleting {}'.format(TmpDir)
+            print('Folder Exists.. Deleting {}'.format(TmpDir))
             shutil.rmtree(TmpDir)
 
         os.mkdir(TmpDir)
diff --git a/scripts/porting_test.sh b/scripts/porting_test.sh
new file mode 100755
index 0000000..23eea19
--- /dev/null
+++ b/scripts/porting_test.sh
@@ -0,0 +1,134 @@
+#!/usr/bin/env bash
+
+# porting_test.sh
+
+
+
+# check the system dependencies for pyprecag are installed.
+# build a temporary virtual environment for python 2 and for
+# python 3 and attempt installation in each.
+
+# set to anything other than "true" to disable verbose testing.
+#VERBOSE="true"
+VERBOSE="false"
+
+if [[ "$1" == "2" ]] ; then
+    TARGET_PYTHONS=("python2.7")
+elif [[ "$1" == "3" ]] ; then
+    TARGET_PYTHONS=("python3.6")
+else
+    TARGET_PYTHONS=("python2.7" "python3.6")
+fi
+
+
+if [[ ${VERBOSE} == "true" ]] ; then
+    VERBOSE_FLAG="-v"
+else
+    VERBOSE_FLAG=""
+fi
+
+
+
+#set -x
+
+function check_requirements {
+    set -e
+    gdal-config --version
+    for P in "${TARGET_PYTHONS[@]}" ; do
+        "${P}" --version
+    done
+    python3 -c 'import virtualenv'
+    echo testing for python3.6-dev
+    python3.6-config --configdir
+    echo OK
+    echo checking we are in a directory with setup.py
+    test -r ./setup.py
+    echo OK
+    set +e
+
+}
+
+function install_gdal {
+    _pip=$1
+    set -x
+    "${_pip}" install GDAL=="$(gdal-config --version | sed -e 's/$//')" --global-option=build_ext --global-option=-I/usr/include/gdal
+    set +x
+    echo "GDAL installation completed."
+    "${_pip}" freeze
+
+}
+
+function install_pvtools {
+    _pip=$1
+    "${_pip}" install -e .[dev,test] || exit 1
+}
+
+
+function run_tests {
+    _python_version=$1
+    venv_directory=".venv/${_python_version}"
+    _python="${venv_directory}/bin/python"
+    _pip="${venv_directory}/bin/pip"
+    "${_pip}" freeze
+    "${_python}" --version
+    (
+        set -x
+        timeout 600s \
+            time "${_python}" -m unittest discover -s pyprecag/tests \
+            -p "*test*.py" "${VERBOSE_FLAG}"
+    )
+}
+
+function run_single_tests {
+    # we have a test (or possibly setup) that is hanging. This is to try and identify which one
+    _python=$1
+    for _test in $(find pyprecag/tests/ -name 'test_*.py') ; do
+        echo testing "${_test}"
+        "${_python}" -m unittest "${_test}"
+    done
+
+}
+
+
+function build_virtualenv {
+  # build a virtualenv (and install dependencies) for a given python, if it doesn't already exist.
+  mkdir -p ".venv"
+  PYTHON=$1
+  venv_directory=".venv/${PYTHON}"
+  _pip="${venv_directory}/bin/pip"
+  echo "testing python virtualenv directory ${venv_directory}"
+  if  test ! -r "${venv_directory}/bin/activate"; then
+      echo virtual environment not found, creating
+      python3 -m virtualenv -p "${PYTHON}" ".venv/${PYTHON}"
+      "${_pip}" install --upgrade pip
+      install_gdal "${_pip}" >/dev/null
+      install_pvtools "${_pip}" >/dev/null
+  else
+      echo "virtual env found; not installing anything."
+  fi
+  ${_pip} freeze
+}
+
+
+check_requirements
+
+TEMPDIR=$(mktemp -d)
+
+for PYTHON in "${TARGET_PYTHONS[@]}" ; do
+  #python3 -m virtualenv -p "${PYTHON}" "${TEMPDIR}/${PYTHON}"
+  _pip="${TEMPDIR}/${PYTHON}/bin/pip"
+  # ls -l "${_pip}"
+  #"${_pip}" --version
+  #install_gdal "${_pip}" >/dev/null
+  #install_pvtools "${_pip}" >/dev/null
+
+
+  _python="${TEMPDIR}/${PYTHON}/bin/python"
+
+  build_virtualenv "${PYTHON}"
+  #run_tests "${_python}"
+  run_tests "${PYTHON}"
+
+done
+
+#find "${TEMPDIR}" -maxdepth 2