From f85ba527f8433447c511cf281ca954e463446c3c Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Tue, 7 Nov 2023 14:03:31 +1100
Subject: [PATCH 1/9] compression with zst and significant digits quantization

---
 gplately/grids.py | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/gplately/grids.py b/gplately/grids.py
index adf34e4a..caa25d3b 100644
--- a/gplately/grids.py
+++ b/gplately/grids.py
@@ -206,6 +206,10 @@ def find_label(keys, labels):
         cdf_lon  = cdf[key_lon][:]
         cdf_lat  = cdf[key_lat][:]
 
+        fill_value = cdf[key_z].missing_value
+
+        cdf_grid[np.isclose(cdf_grid, fill_value, rtol=0.1)] = np.nan
+
     if realign:
         # realign longitudes to -180/180 dateline
         cdf_grid_z, cdf_lon, cdf_lat = realign_grid(cdf_grid, cdf_lon, cdf_lat)
@@ -243,7 +247,7 @@ def find_label(keys, labels):
     else:
         return cdf_grid_z
     
-def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90]):
+def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90], significant_digits=None, fill_value=np.nan):
     """ Write geological data contained in a `grid` to a netCDF4 grid with a specified `filename`.
 
     Notes
@@ -288,8 +292,8 @@ def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90]):
         cdf.title = "Grid produced by gplately"
         cdf.createDimension('lon', lon_grid.size)
         cdf.createDimension('lat', lat_grid.size)
-        cdf_lon = cdf.createVariable('lon', lon_grid.dtype, ('lon',), zlib=True)
-        cdf_lat = cdf.createVariable('lat', lat_grid.dtype, ('lat',), zlib=True)
+        cdf_lon = cdf.createVariable('lon', lon_grid.dtype, ('lon',), compression='zstd', complevel=9)
+        cdf_lat = cdf.createVariable('lat', lat_grid.dtype, ('lat',), compression='zstd', complevel=9)
         cdf_lon[:] = lon_grid
         cdf_lat[:] = lat_grid
 
@@ -301,10 +305,18 @@ def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90]):
         cdf_lat.standard_name = 'lat'
         cdf_lat.actual_range = [lat_grid[0], lat_grid[-1]]
 
-        cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), zlib=True)
+        if significant_digits:
+            cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), compression='zstd', complevel=9,
+                                          significant_digits=int(significant_digits),
+                                          quantize_mode='GranularBitRound')
+
+        else:
+            cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), compression='zstd', complevel=9)
+
         # netCDF4 uses the missing_value attribute as the default _FillValue
         # without this, _FillValue defaults to 9.969209968386869e+36
-        cdf_data.missing_value = np.nan
+        cdf_data.missing_value = fill_value
+
         cdf_data.standard_name = 'z'
         #Ensure pygmt registers min and max z values properly
         cdf_data.actual_range = [np.nanmin(grid), np.nanmax(grid)]

From d2e8339ed0d58567e50fb3f424e0b235b2a5b793 Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Tue, 7 Nov 2023 14:29:58 +1100
Subject: [PATCH 2/9] read netcdf files even when there is no missing_values
 attribute

---
 gplately/grids.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/gplately/grids.py b/gplately/grids.py
index caa25d3b..b2d857b8 100644
--- a/gplately/grids.py
+++ b/gplately/grids.py
@@ -206,9 +206,9 @@ def find_label(keys, labels):
         cdf_lon  = cdf[key_lon][:]
         cdf_lat  = cdf[key_lat][:]
 
-        fill_value = cdf[key_z].missing_value
-
-        cdf_grid[np.isclose(cdf_grid, fill_value, rtol=0.1)] = np.nan
+        if hasattr(cdf[key_z], 'missing_value'):
+            fill_value = cdf[key_z].missing_value
+            cdf_grid[np.isclose(cdf_grid, fill_value, rtol=0.1)] = np.nan
 
     if realign:
         # realign longitudes to -180/180 dateline

From 8f8102a0a3b40ac9c469633f771aa5106fde6f56 Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Fri, 13 Sep 2024 21:39:53 +1000
Subject: [PATCH 3/9] update write_netcdf_grid

---
 gplately/grids.py | 41 +++++++++++++++++++++++++++++++----------
 1 file changed, 31 insertions(+), 10 deletions(-)

diff --git a/gplately/grids.py b/gplately/grids.py
index b2d857b8..fad72a54 100644
--- a/gplately/grids.py
+++ b/gplately/grids.py
@@ -247,7 +247,7 @@ def find_label(keys, labels):
     else:
         return cdf_grid_z
     
-def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90], significant_digits=None, fill_value=np.nan):
+def write_netcdf_grid(filename, grid, extent="global", significant_digits=None, fill_value=np.nan, **kwargs):
     """ Write geological data contained in a `grid` to a netCDF4 grid with a specified `filename`.
 
     Notes
@@ -272,7 +272,7 @@ def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90], significant_digi
         An ndarray grid containing data to be written into a `netCDF` (.nc) file. Note: Rows correspond to 
         the data's latitudes, while the columns correspond to the data's longitudes.
 
-    extent : 1D numpy array, default=[-180,180,-90,90]
+    extent : list, default=[-180,180,-90,90]
         Four elements that specify the [min lon, max lon, min lat, max lat] to constrain the lat and lon 
         variables of the netCDF grid to. If no extents are supplied, full global extent `[-180, 180, -90, 90]` 
         is assumed. 
@@ -282,18 +282,25 @@ def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90], significant_digi
     A netCDF grid will be saved to the path specified in `filename`.
     """
     import netCDF4
+
+    if extent == 'global':
+        extent = [-180, 180, -90, 90]
+    else:
+        assert len(extent) == 4, "specify the [min lon, max lon, min lat, max lat]"
     
     nrows, ncols = np.shape(grid)
     
     lon_grid = np.linspace(extent[0], extent[1], ncols)
     lat_grid = np.linspace(extent[2], extent[3], nrows)
+
+    data_kwds = {'compression': 'zlib', 'complevel': 9}
     
     with netCDF4.Dataset(filename, 'w', driver=None) as cdf:
         cdf.title = "Grid produced by gplately"
         cdf.createDimension('lon', lon_grid.size)
         cdf.createDimension('lat', lat_grid.size)
-        cdf_lon = cdf.createVariable('lon', lon_grid.dtype, ('lon',), compression='zstd', complevel=9)
-        cdf_lat = cdf.createVariable('lat', lat_grid.dtype, ('lat',), compression='zstd', complevel=9)
+        cdf_lon = cdf.createVariable('lon', lon_grid.dtype, ('lon',), **data_kwds)
+        cdf_lat = cdf.createVariable('lat', lat_grid.dtype, ('lat',), **data_kwds)
         cdf_lon[:] = lon_grid
         cdf_lat[:] = lat_grid
 
@@ -305,22 +312,36 @@ def write_netcdf_grid(filename, grid, extent=[-180,180,-90,90], significant_digi
         cdf_lat.standard_name = 'lat'
         cdf_lat.actual_range = [lat_grid[0], lat_grid[-1]]
 
+        # create container variable for CRS: lon/lat WGS84 datum
+        crso = cdf.createVariable('crs','i4')
+        crso.long_name = 'Lon/Lat Coords in WGS84'
+        crso.grid_mapping_name='latitude_longitude'
+        crso.longitude_of_prime_meridian = 0.0
+        crso.semi_major_axis = 6378137.0
+        crso.inverse_flattening = 298.257223563
+        crso.spatial_ref = """GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]"""
+
+        # add more keyword arguments for quantizing data
         if significant_digits:
-            cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), compression='zstd', complevel=9,
-                                          significant_digits=int(significant_digits),
-                                          quantize_mode='GranularBitRound')
+            data_kwds['significant_digits'] = int(significant_digits)
+            data_kwds['quantize_mode'] = 'GranularBitRound'
 
-        else:
-            cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), compression='zstd', complevel=9)
+        cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), **data_kwds)
 
         # netCDF4 uses the missing_value attribute as the default _FillValue
         # without this, _FillValue defaults to 9.969209968386869e+36
         cdf_data.missing_value = fill_value
 
         cdf_data.standard_name = 'z'
-        #Ensure pygmt registers min and max z values properly
+
+        cdf_data.add_offset = 0.0
+        cdf_data.grid_mapping = 'crs'
+        cdf_data.set_auto_maskandscale(False)
+        
+        # ensure min and max z values are properly registered
         cdf_data.actual_range = [np.nanmin(grid), np.nanmax(grid)]
 
+        # write data
         cdf_data[:,:] = grid
 
 

From 97fcfa3daa396f5426ddf9e1f42e3dfa230f474c Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Mon, 16 Sep 2024 11:50:48 +1000
Subject: [PATCH 4/9] add default significant digits (2) and descriptive title

---
 gplately/grids.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/gplately/grids.py b/gplately/grids.py
index 30843b96..30cea1e7 100644
--- a/gplately/grids.py
+++ b/gplately/grids.py
@@ -314,6 +314,7 @@ def write_netcdf_grid(filename, grid, extent="global", significant_digits=None,
     A netCDF grid will be saved to the path specified in `filename`.
     """
     import netCDF4
+    from gplately import __version__ as _version
 
     if extent == 'global':
         extent = [-180, 180, -90, 90]
@@ -328,7 +329,7 @@ def write_netcdf_grid(filename, grid, extent="global", significant_digits=None,
     data_kwds = {'compression': 'zlib', 'complevel': 9}
     
     with netCDF4.Dataset(filename, 'w', driver=None) as cdf:
-        cdf.title = "Grid produced by gplately"
+        cdf.title = "Grid produced by gplately " + str(_version)
         cdf.createDimension('lon', lon_grid.size)
         cdf.createDimension('lat', lat_grid.size)
         cdf_lon = cdf.createVariable('lon', lon_grid.dtype, ('lon',), **data_kwds)
@@ -355,7 +356,8 @@ def write_netcdf_grid(filename, grid, extent="global", significant_digits=None,
 
         # add more keyword arguments for quantizing data
         if significant_digits:
-            data_kwds['significant_digits'] = int(significant_digits)
+            # significant_digits needs to be >= 2 so that NaNs are preserved
+            data_kwds['significant_digits'] = max(2, int(significant_digits))
             data_kwds['quantize_mode'] = 'GranularBitRound'
 
         cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), **data_kwds)

From 78485e3e545696b47d9db85713c074970b3e41a7 Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Mon, 16 Sep 2024 15:50:17 +1000
Subject: [PATCH 5/9] add support for boolean grids

---
 gplately/grids.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/gplately/grids.py b/gplately/grids.py
index 30cea1e7..631dcf4d 100644
--- a/gplately/grids.py
+++ b/gplately/grids.py
@@ -360,11 +360,18 @@ def write_netcdf_grid(filename, grid, extent="global", significant_digits=None,
             data_kwds['significant_digits'] = max(2, int(significant_digits))
             data_kwds['quantize_mode'] = 'GranularBitRound'
 
+        # boolean arrays need to be converted to integers
+        # no such thing as a mask on a boolean array
+        if grid.dtype is np.dtype(bool):
+            grid = grid.astype('i1')
+            fill_value = None
+
         cdf_data = cdf.createVariable('z', grid.dtype, ('lat','lon'), **data_kwds)
 
         # netCDF4 uses the missing_value attribute as the default _FillValue
         # without this, _FillValue defaults to 9.969209968386869e+36
-        cdf_data.missing_value = fill_value
+        if fill_value is not None:
+            cdf_data.missing_value = fill_value
 
         cdf_data.standard_name = 'z'
 

From 91008d10b9116c15f4783ac6b9697a011c0b6f70 Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Mon, 16 Sep 2024 15:50:45 +1000
Subject: [PATCH 6/9] update age gridding workflow to spit out optimised netcdf
 grids

---
 gplately/oceans.py | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/gplately/oceans.py b/gplately/oceans.py
index d1dec4e7..3e3c81cd 100644
--- a/gplately/oceans.py
+++ b/gplately/oceans.py
@@ -781,8 +781,9 @@ def _create_continental_mask(self, time_array):
 
             grids.write_netcdf_grid(
                 self.continent_mask_filepath.format(time),
-                final_grid,
+                final_grid.astype('i1'),
                 extent=[-180, 180, -90, 90],
+                fill_value=None,
             )
             logger.info(f"Finished building a continental mask at {time} Ma!")
 
@@ -809,8 +810,9 @@ def _build_continental_mask(self, time: float, overwrite=False):
 
         grids.write_netcdf_grid(
             self.continent_mask_filepath.format(time),
-            final_grid,
+            final_grid.astype('i1'),
             extent=[-180, 180, -90, 90],
+            fill_value=None,
         )
         logger.info(f"Finished building a continental mask at {time} Ma!")
 
@@ -1490,7 +1492,11 @@ def _save_netcdf_file(
     grid_output = os.path.join(output_dir, grid_basename)
 
     if unmasked:
-        grids.write_netcdf_grid(grid_output_unmasked, Z, extent=extent)
+        grids.write_netcdf_grid(
+            grid_output_unmasked, 
+            Z, 
+            extent=extent,
+            significant_digits=2)
 
     # Identify regions in the grid in the continental mask
     cont_mask = grids.Raster(data=continent_mask_filename.format(time))
@@ -1512,6 +1518,7 @@ def _save_netcdf_file(
         grid_output,
         Z,
         extent=extent,
+        significant_digits=2,
     )
     logger.info(f"Save {name} netCDF grid at {time:0.2f} Ma completed!")
 
@@ -1578,7 +1585,10 @@ def _lat_lon_z_to_netCDF_time(
     grid_output = os.path.join(output_dir, grid_basename)
 
     if unmasked:
-        grids.write_netcdf_grid(grid_output_unmasked, Z, extent=extent)
+        grids.write_netcdf_grid(grid_output_unmasked, 
+            Z, 
+            extent=extent,
+            significant_digits=2)
 
     # Identify regions in the grid in the continental mask
     cont_mask = grids.Raster(data=continent_mask_filename.format(time))
@@ -1603,6 +1613,7 @@ def _lat_lon_z_to_netCDF_time(
         grid_output,
         Z,
         extent=extent,
+        significant_digits=2,
     )
     logger.info(f"{zval_name} netCDF grids for {time:0.2f} Ma complete!")
 
@@ -1717,8 +1728,9 @@ def continent_contouring_buffer_and_gap_distance_radians(time, contoured_contine
     )
     grids.write_netcdf_grid(
         continent_mask_filepath.format(time),
-        continent_mask.astype("float"),
+        continent_mask.astype("i1"),
         extent=[-180, 180, -90, 90],
+        fill_value=None,
     )
     logger.info(
         f"Finished building a continental mask at {time} Ma! (continent_contouring)"

From d02915be264d56eea6519ce16f8b9235a9c46f70 Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Tue, 17 Sep 2024 10:28:35 +1000
Subject: [PATCH 7/9] add resize parameter for reading netcdf grids

---
 gplately/grids.py | 100 +++++++++++++++++++++++++++++++++++-----------
 1 file changed, 77 insertions(+), 23 deletions(-)

diff --git a/gplately/grids.py b/gplately/grids.py
index 631dcf4d..a8791821 100644
--- a/gplately/grids.py
+++ b/gplately/grids.py
@@ -142,7 +142,7 @@ def realign_grid(array, lons, lats):
     return array, lons, lats
 
 
-def read_netcdf_grid(filename, return_grids=False, realign=False, resample=None):
+def read_netcdf_grid(filename, return_grids=False, realign=False, resample=None, resize=None):
     """Read a `netCDF` (.nc) grid from a given `filename` and return its data as a
     `MaskedArray`.
 
@@ -174,6 +174,10 @@ def read_netcdf_grid(filename, return_grids=False, realign=False, resample=None)
         If passed as `resample = (spacingX, spacingY)`, the given `netCDF` grid is resampled
         with these x and y resolutions.
 
+    resize : tuple, optional, default=None
+        If passed as `resample = (resX, resY)`, the given `netCDF` grid is resized
+        to the number of columns (resX) and rows (resY).
+
     Returns
     -------
     grid_z : MaskedArray
@@ -237,10 +241,18 @@ def find_label(keys, labels):
         cdf_lon = cdf[key_lon][:]
         cdf_lat = cdf[key_lat][:]
 
-        if hasattr(cdf[key_z], 'missing_value'):
+        # fill missing values
+        if hasattr(cdf[key_z], 'missing_value') and np.issubdtype(cdf_grid.dtype, np.floating):
             fill_value = cdf[key_z].missing_value
             cdf_grid[np.isclose(cdf_grid, fill_value, rtol=0.1)] = np.nan
 
+        # convert to boolean array
+        if np.issubdtype(cdf_grid.dtype, np.integer):
+            unique_grid = np.unique(cdf_grid)
+            if len(unique_grid) == 2:
+                if (unique_grid == [0,1]).all():
+                    cdf_grid = cdf_grid.astype(bool)
+
     if realign:
         # realign longitudes to -180/180 dateline
         cdf_grid_z, cdf_lon, cdf_lat = realign_grid(cdf_grid, cdf_lon, cdf_lat)
@@ -250,25 +262,58 @@ def find_label(keys, labels):
     # resample
     if resample is not None:
         spacingX, spacingY = resample
-        lon_grid = np.arange(cdf_lon.min(), cdf_lon.max() + spacingX, spacingX)
-        lat_grid = np.arange(cdf_lat.min(), cdf_lat.max() + spacingY, spacingY)
-        lonq, latq = np.meshgrid(lon_grid, lat_grid)
-        original_extent = (
-            cdf_lon[0],
-            cdf_lon[-1],
-            cdf_lat[0],
-            cdf_lat[-1],
-        )
-        cdf_grid_z = sample_grid(
-            lonq,
-            latq,
-            cdf_grid_z,
-            method="nearest",
-            extent=original_extent,
-            return_indices=False,
-        )
-        cdf_lon = lon_grid
-        cdf_lat = lat_grid
+
+        # don't resample if already the same resolution
+        dX = np.diff(cdf_lon).mean()
+        dY = np.diff(cdf_lat).mean()
+
+        if spacingX != dX or spacingY != dY:
+            lon_grid = np.arange(cdf_lon.min(), cdf_lon.max() + spacingX, spacingX)
+            lat_grid = np.arange(cdf_lat.min(), cdf_lat.max() + spacingY, spacingY)
+            lonq, latq = np.meshgrid(lon_grid, lat_grid)
+            original_extent = (
+                cdf_lon[0],
+                cdf_lon[-1],
+                cdf_lat[0],
+                cdf_lat[-1],
+            )
+            cdf_grid_z = sample_grid(
+                lonq,
+                latq,
+                cdf_grid_z,
+                method="nearest",
+                extent=original_extent,
+                return_indices=False,
+            )
+            cdf_lon = lon_grid
+            cdf_lat = lat_grid
+
+    # resize
+    if resize is not None:
+        resX, resY = resize
+
+        # don't resize if already the same shape
+        if resX != cdf_grid_z.shape[1] or resY != cdf_grid_z.shape[0]:
+            original_extent = (
+                cdf_lon[0],
+                cdf_lon[-1],
+                cdf_lat[0],
+                cdf_lat[-1],
+            )
+            lon_grid = np.linspace(original_extent[0], original_extent[1], resX)
+            lat_grid = np.linspace(original_extent[2], original_extent[3], resY)
+            lonq, latq = np.meshgrid(lon_grid, lat_grid)
+
+            cdf_grid_z = sample_grid(
+                lonq,
+                latq,
+                cdf_grid_z,
+                method="nearest",
+                extent=original_extent,
+                return_indices=False,
+            )
+            cdf_lon = lon_grid
+            cdf_lat = lat_grid
 
     # Fix grids with 9e36 as the fill value for nan.
     # cdf_grid_z.fill_value = float('nan')
@@ -1501,6 +1546,7 @@ def __init__(
         extent="global",
         realign=False,
         resample=None,
+        resize=None,
         time=0.0,
         origin=None,
         **kwargs,
@@ -1530,6 +1576,10 @@ def __init__(
             Optionally resample grid, pass spacing in X and Y direction as a
             2-tuple e.g. resample=(spacingX, spacingY).
 
+        resize : 2-tuple, optional
+            Optionally resample grid to X-columns, Y-rows as a
+            2-tuple e.g. resample=(resX, resY).
+
         time : float, default: 0.0
             The time step represented by the raster data. Used for raster
             reconstruction.
@@ -1600,6 +1650,7 @@ def __init__(
                 return_grids=True,
                 realign=realign,
                 resample=resample,
+                resize=resize,
             )
             self._lons = lons
             self._lats = lats
@@ -1621,6 +1672,9 @@ def __init__(
         if (not isinstance(data, str)) and (resample is not None):
             self.resample(*resample, inplace=True)
 
+        if (not isinstance(data, str)) and (resize is not None):
+            self.resize(*resize, inplace=True)
+
     @property
     def time(self):
         """The time step represented by the raster data."""
@@ -1963,10 +2017,10 @@ def fill_NaNs(self, inplace=False, return_array=False):
         else:
             return Raster(data, self.plate_reconstruction, self.extent, self.time)
 
-    def save_to_netcdf4(self, filename):
+    def save_to_netcdf4(self, filename, significant_digits=None, fill_value=np.nan):
         """Saves the grid attributed to the `Raster` object to the given `filename` (including
         the ".nc" extension) in netCDF4 format."""
-        write_netcdf_grid(str(filename), self.data, self.extent)
+        write_netcdf_grid(str(filename), self.data, self.extent, significant_digits, fill_value)
 
     def reconstruct(
         self,

From b635a23964ac0b2a499b7048c6851683c082b809 Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Tue, 17 Sep 2024 10:29:04 +1000
Subject: [PATCH 8/9] refactor writing masked and unmasked grids

---
 gplately/oceans.py | 38 ++++++++++++++------------------------
 1 file changed, 14 insertions(+), 24 deletions(-)

diff --git a/gplately/oceans.py b/gplately/oceans.py
index 3e3c81cd..f3aaaf5f 100644
--- a/gplately/oceans.py
+++ b/gplately/oceans.py
@@ -1480,6 +1480,7 @@ def _save_netcdf_file(
 
     # Interpolate lons, lats and zvals over a regular grid using nearest neighbour interpolation
     Z = tools.griddata_sphere((lons, lats), zdata, (X, Y), method="nearest")
+    Z = Z.astype('f4')
 
     unmasked_basename = f"{name}_grid_unmasked_{time:0.2f}_Ma.nc"
     grid_basename = f"{name}_grid_{time:0.2f}_Ma.nc"
@@ -1496,29 +1497,24 @@ def _save_netcdf_file(
             grid_output_unmasked, 
             Z, 
             extent=extent,
-            significant_digits=2)
+            significant_digits=2,
+            fill_value=None)
 
     # Identify regions in the grid in the continental mask
-    cont_mask = grids.Raster(data=continent_mask_filename.format(time))
-
     # We need the continental mask to match the number of nodes
     # in the uniform grid defined above. This is important if we
     # pass our own continental mask to SeafloorGrid
-    if cont_mask.shape[1] != resX:
-        cont_mask.resize(resX, resY, inplace=True)
+    cont_mask = grids.read_netcdf_grid(continent_mask_filename.format(time), resize=(resX,resY))
 
     # Use the continental mask
-    Z = np.ma.array(
-        grids.Raster(data=Z.astype("float")).data.data,
-        mask=cont_mask.data.data,
-        fill_value=np.nan,
-    )
+    Z[cont_mask] = np.nan
 
     grids.write_netcdf_grid(
         grid_output,
         Z,
         extent=extent,
         significant_digits=2,
+        fill_value=np.nan,
     )
     logger.info(f"Save {name} netCDF grid at {time:0.2f} Ma completed!")
 
@@ -1573,6 +1569,7 @@ def _lat_lon_z_to_netCDF_time(
 
     # Interpolate lons, lats and zvals over a regular grid using nearest neighbour interpolation
     Z = tools.griddata_sphere((lons, lats), zdata, (X, Y), method="nearest")
+    Z = Z.astype('f4') # float32 precision
 
     unmasked_basename = f"{zval_name}_grid_unmasked_{time:0.2f}Ma.nc"
     grid_basename = f"{zval_name}_grid_{time:0.2f}Ma.nc"
@@ -1585,35 +1582,28 @@ def _lat_lon_z_to_netCDF_time(
     grid_output = os.path.join(output_dir, grid_basename)
 
     if unmasked:
-        grids.write_netcdf_grid(grid_output_unmasked, 
+        grids.write_netcdf_grid(
+            grid_output_unmasked, 
             Z, 
             extent=extent,
-            significant_digits=2)
+            significant_digits=2,
+            fill_value=None)
 
     # Identify regions in the grid in the continental mask
-    cont_mask = grids.Raster(data=continent_mask_filename.format(time))
-
     # We need the continental mask to match the number of nodes
     # in the uniform grid defined above. This is important if we
     # pass our own continental mask to SeafloorGrid
-    if cont_mask.shape[1] != resX:
-        cont_mask.resize(resX, resY, inplace=True)
+    cont_mask = grids.read_netcdf_grid(continent_mask_filename.format(time), resize=(resX,resY))
 
     # Use the continental mask
-    Z = np.ma.array(
-        grids.Raster(data=Z.astype("float")).data.data,
-        mask=cont_mask.data.data,
-        fill_value=np.nan,
-    )
-
-    # grd = cont_mask.interpolate(X, Y) > 0.5
-    # Z[grd] = np.nan
+    Z[cont_mask] = np.nan
 
     grids.write_netcdf_grid(
         grid_output,
         Z,
         extent=extent,
         significant_digits=2,
+        fill_value=np.nan,
     )
     logger.info(f"{zval_name} netCDF grids for {time:0.2f} Ma complete!")
 

From b719ab3ac47dfdad264a51d3e4b0029e25e0f48a Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Tue, 17 Sep 2024 11:17:43 +1000
Subject: [PATCH 9/9] moving code from reconstruction.py to
 seafloor_grid_utils.py

---
 gplately/__init__.py                  |  10 -
 gplately/oceans.py                    |   6 +-
 gplately/reconstruction.py            | 785 -------------------------
 gplately/utils/__init__.py            |   8 +
 gplately/utils/seafloor_grid_utils.py | 787 +++++++++++++++++++++++++-
 5 files changed, 796 insertions(+), 800 deletions(-)

diff --git a/gplately/__init__.py b/gplately/__init__.py
index 98255dfe..fb83a9ca 100644
--- a/gplately/__init__.py
+++ b/gplately/__init__.py
@@ -232,7 +232,6 @@
     pygplates,
     reconstruction,
 )
-from .data import DataCollection
 from .download import DataServer
 from .grids import Raster
 from .oceans import SeafloorGrid
@@ -240,9 +239,6 @@
 from .reconstruction import (
     PlateReconstruction,
     Points,
-    _ContinentCollision,
-    _DefaultCollision,
-    _ReconstructByTopologies,
 )
 from .tools import EARTH_RADIUS
 from .utils import io_utils
@@ -250,9 +246,6 @@
 
 __pdoc__ = {
     "data": False,
-    "_DefaultCollision": False,
-    "_ContinentCollision": False,
-    "_ReconstructByTopologies": False,
     "examples": False,
     "notebooks": False,
     "commands": False,
@@ -283,9 +276,6 @@
     "Points",
     "Raster",
     "SeafloorGrid",
-    "_ContinentCollision",
-    "_DefaultCollision",
-    "_ReconstructByTopologies",
     # Functions
     "get_geometries",
     "get_valid_geometries",
diff --git a/gplately/oceans.py b/gplately/oceans.py
index f3aaaf5f..9de9b8c8 100644
--- a/gplately/oceans.py
+++ b/gplately/oceans.py
@@ -1147,7 +1147,7 @@ def reconstruct_by_topologies(self):
         point_id = range(len(active_points))
 
         # Specify the default collision detection region as subduction zones
-        default_collision = reconstruction._DefaultCollision(
+        default_collision = seafloor_grid_utils.DefaultCollision(
             feature_specific_collision_parameters=[
                 (
                     pygplates.FeatureType.gpml_subduction_zone,
@@ -1156,7 +1156,7 @@ def reconstruct_by_topologies(self):
             ]
         )
         # In addition to the default subduction detection, also detect continental collisions
-        collision_spec = reconstruction._ContinentCollision(
+        collision_spec = seafloor_grid_utils.ContinentCollision(
             # This filename string should not have a time formatted into it - this is
             # taken care of later.
             self.continent_mask_filepath,
@@ -1165,7 +1165,7 @@ def reconstruct_by_topologies(self):
         )
 
         # Call the reconstruct by topologies object
-        topology_reconstruction = reconstruction._ReconstructByTopologies(
+        topology_reconstruction = seafloor_grid_utils.ReconstructByTopologies(
             self.rotation_model,
             self.topology_features,
             self._max_time,
diff --git a/gplately/reconstruction.py b/gplately/reconstruction.py
index ad383576..262b2cd9 100644
--- a/gplately/reconstruction.py
+++ b/gplately/reconstruction.py
@@ -20,8 +20,6 @@
 working with point data, and calculating plate velocities at specific geological times. 
 """
 
-import math
-import os
 import warnings
 
 import numpy as np
@@ -2009,786 +2007,3 @@ def rotate_reference_frames(
                 time=reconstruction_time,
                 plate_id=self.plate_id,
             )
-
-
-# FROM RECONSTRUCT_BY_TOPOLOGIES.PY
-class _DefaultCollision(object):
-    """
-    Default collision detection function class (the function is the '__call__' method).
-    """
-
-    DEFAULT_GLOBAL_COLLISION_PARAMETERS = (7.0, 10.0)
-    """
-    Default collision parameters for all feature types.
-
-    This is a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my):
-    Here we default to the same constants used internally in GPlates 2.0 (ie, 7.0 and 10.0).
-    """
-
-    def __init__(
-        self,
-        global_collision_parameters=DEFAULT_GLOBAL_COLLISION_PARAMETERS,
-        feature_specific_collision_parameters=None,
-    ):
-        """
-        global_collision_parameters: The collision parameters to use for any feature type not specified in 'feature_specific_collision_parameters'.
-                                     Should be a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my).
-                                     The first threshold parameter means:
-                                        A point that transitions from one plate to another can disappear if the change in velocity exceeds this threshold.
-                                     The second threshold parameter means:
-                                        Only those transitioning points exceeding the threshold velocity delta and that are close enough to a plate boundary can disappear.
-                                        The distance is proportional to the relative velocity (change in velocity), plus a constant offset based on the threshold distance to boundary
-                                        to account for plate boundaries that change shape significantly from one time step to the next
-                                        (note that some boundaries are meant to do this and others are a result of digitisation).
-                                        The actual distance threshold used is (threshold_distance_to_boundary + relative_velocity) * time_interval
-                                     Defaults to parameters used in GPlates 2.0, if not specified.
-
-        feature_specific_collision_parameters: Optional sequence of collision parameters specific to feature types.
-                                               If specified then should be a sequence of 2-tuples, with each 2-tuple specifying (feature_type, collision_parameters).
-                                               And where each 'collision_parameters' is a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my).
-                                                   See 'global_collision_parameters' for details on these thresholds.
-                                               Any feature type not specified here defaults to using 'global_collision_parameters'.
-        """
-
-        # Convert list of (feature_type, collision_parameters) tuples to a dictionary.
-        if feature_specific_collision_parameters:
-            self.feature_specific_collision_parameters = dict(
-                feature_specific_collision_parameters
-            )
-        else:
-            self.feature_specific_collision_parameters = dict()
-        # Fallback for any feature type not specified in the optional feature-specific list.
-        self.global_collision_parameters = global_collision_parameters
-
-        # Used to improve performance by caching velocity stage rotations in a dict (for a specific reconstruction time).
-        self.velocity_stage_rotation_dict = {}
-        self.velocity_stage_rotation_time = None
-
-    def __call__(
-        self,
-        rotation_model,
-        time,
-        reconstruction_time_interval,
-        prev_point,
-        curr_point,
-        prev_topology_plate_id,
-        prev_resolved_plate_boundary,
-        curr_topology_plate_id,
-        curr_resolved_plate_boundary,
-    ):
-        """
-        Returns True if a collision was detected.
-
-        If transitioning from a rigid plate to another rigid plate with a different plate ID then
-        calculate the difference in velocities and continue testing as follows
-        (otherwise, if there's no transition, then the point is still active)...
-
-        If the velocity difference is below a threshold then we assume the previous plate was split,
-        or two plates joined. In this case the point has not subducted (forward in time) or
-        been consumed by a mid-ocean (backward in time) and hence is still active.
-
-        If the velocity difference is large enough then we see if the distance of the *previous* position
-        to the polygon boundary (of rigid plate containing it) exceeds a threshold.
-        If the distance exceeds the threshold then the point is far enough away from the boundary that it
-        cannot be subducted or consumed by it and hence the point is still active.
-        However if the point is close enough then we assume the point was subducted/consumed
-        (remember that the point switched plate IDs).
-        Also note that the threshold distance increases according to the velocity difference to account for fast
-        moving points (that would otherwise tunnel through the boundary and accrete onto the other plate).
-        The reason for testing the distance from the *previous* point, and not from the *current* point, is:
-
-          (i)  A topological boundary may *appear* near the current point (such as a plate split at the current time)
-               and we don't want that split to consume the current point regardless of the velocity difference.
-               It won't get consumed because the *previous* point was not near a boundary (because before split happened).
-               If the velocity difference is large enough then it might cause the current point to transition to the
-               adjacent split plate in the *next* time step (and that's when it should get consumed, not in the current time step).
-               An example of this is a mid-ocean ridge suddenly appearing (going forward in time).
-
-          (ii) A topological boundary may *disappear* near the current point (such as a plate merge at the current time)
-               and we want that merge to consume the current point if the velocity difference is large enough.
-               In this case the *previous* point is near a boundary (because before plate merged) and hence can be
-               consumed (provided velocity difference is large enough). And since the boundary existed in the previous
-               time step, it will affect position of the current point (and whether it gets consumed or not).
-               An example of this is a mid-ocean ridge suddenly disappearing (going backward in time).
-
-        ...note that items (i) and (ii) above apply both going forward and backward in time.
-        """
-
-        # See if a collision occurred.
-        if (
-            curr_topology_plate_id != prev_topology_plate_id
-            and prev_topology_plate_id is not None
-            and curr_topology_plate_id is not None
-        ):
-            #
-            # Speed up by caching velocity stage rotations in a dict.
-            #
-            if time != self.velocity_stage_rotation_time:
-                # We've just switched to a new time so clear the cache.
-                #
-                # We only cache stage rotations for a specific time.
-                # We only really need to cache different plate IDs at the same 'time', so this avoids caching for all times
-                # (which would also require including 'time' in the key) and using memory unnecessarily.
-                self.velocity_stage_rotation_dict.clear()
-                self.velocity_stage_rotation_time = time
-            prev_location_velocity_stage_rotation = (
-                self.velocity_stage_rotation_dict.get(prev_topology_plate_id)
-            )
-            if not prev_location_velocity_stage_rotation:
-                prev_location_velocity_stage_rotation = rotation_model.get_rotation(
-                    time + 1, prev_topology_plate_id, time
-                )
-                self.velocity_stage_rotation_dict[prev_topology_plate_id] = (
-                    prev_location_velocity_stage_rotation
-                )
-            curr_location_velocity_stage_rotation = (
-                self.velocity_stage_rotation_dict.get(curr_topology_plate_id)
-            )
-            if not curr_location_velocity_stage_rotation:
-                curr_location_velocity_stage_rotation = rotation_model.get_rotation(
-                    time + 1, curr_topology_plate_id, time
-                )
-                self.velocity_stage_rotation_dict[curr_topology_plate_id] = (
-                    curr_location_velocity_stage_rotation
-                )
-
-            # Note that even though the current point is not inside the previous boundary (because different plate ID), we can still
-            # calculate a velocity using its plate ID (because we really should use the same point in our velocity comparison).
-            prev_location_velocity = pygplates.calculate_velocities(
-                (curr_point,),
-                prev_location_velocity_stage_rotation,
-                1,
-                pygplates.VelocityUnits.kms_per_my,
-            )[0]
-            curr_location_velocity = pygplates.calculate_velocities(
-                (curr_point,),
-                curr_location_velocity_stage_rotation,
-                1,
-                pygplates.VelocityUnits.kms_per_my,
-            )[0]
-
-            delta_velocity = curr_location_velocity - prev_location_velocity
-            delta_velocity_magnitude = delta_velocity.get_magnitude()
-
-            # If we have feature-specific collision parameters then iterate over the boundary sub-segments of the *previous* topological boundary
-            # and test proximity to each sub-segment individually (with sub-segment feature type specific collision parameters).
-            # Otherwise just test proximity to the entire boundary polygon using the global collision parameters.
-            if self.feature_specific_collision_parameters:
-                for (
-                    prev_boundary_sub_segment
-                ) in prev_resolved_plate_boundary.get_boundary_sub_segments():
-                    # Use feature-specific collision parameters if found (falling back to global collision parameters).
-                    (
-                        threshold_velocity_delta,
-                        threshold_distance_to_boundary_per_my,
-                    ) = self.feature_specific_collision_parameters.get(
-                        prev_boundary_sub_segment.get_feature().get_feature_type(),
-                        # Default to global collision parameters if no collision parameters specified for sub-segment's feature type...
-                        self.global_collision_parameters,
-                    )
-
-                    # Since each feature type could use different collision parameters we must use the current boundary sub-segment instead of the boundary polygon.
-                    if self._detect_collision_using_collision_parameters(
-                        reconstruction_time_interval,
-                        delta_velocity_magnitude,
-                        prev_point,
-                        prev_boundary_sub_segment.get_resolved_geometry(),
-                        threshold_velocity_delta,
-                        threshold_distance_to_boundary_per_my,
-                    ):
-                        # Detected a collision.
-                        return True
-            else:
-                # No feature-specific collision parameters so use global fallback.
-                (
-                    threshold_velocity_delta,
-                    threshold_distance_to_boundary_per_my,
-                ) = self.global_collision_parameters
-
-                # Since all feature types use the same collision parameters we can use the boundary polygon instead of iterating over its sub-segments.
-                if self._detect_collision_using_collision_parameters(
-                    reconstruction_time_interval,
-                    delta_velocity_magnitude,
-                    prev_point,
-                    prev_resolved_plate_boundary.get_resolved_boundary(),
-                    threshold_velocity_delta,
-                    threshold_distance_to_boundary_per_my,
-                ):
-                    # Detected a collision.
-                    return True
-
-        return False
-
-    def _detect_collision_using_collision_parameters(
-        self,
-        reconstruction_time_interval,
-        delta_velocity_magnitude,
-        prev_point,
-        prev_boundary_geometry,
-        threshold_velocity_delta,
-        threshold_distance_to_boundary_per_my,
-    ):
-        if delta_velocity_magnitude > threshold_velocity_delta:
-            # Add the minimum distance threshold to the delta velocity threshold.
-            # The delta velocity threshold only allows those points that are close enough to the boundary to reach
-            # it given their current relative velocity.
-            # The minimum distance threshold accounts for sudden changes in the shape of a plate boundary
-            # which are no supposed to represent a new or shifted boundary but are just a result of the topology
-            # builder/user digitising a new boundary line that differs noticeably from that of the previous time period.
-            distance_threshold_radians = (
-                (threshold_distance_to_boundary_per_my + delta_velocity_magnitude)
-                * reconstruction_time_interval
-                / pygplates.Earth.equatorial_radius_in_kms
-            )
-            distance_threshold_radians = min(distance_threshold_radians, math.pi)
-            distance_threshold_radians = max(distance_threshold_radians, 0.0)
-
-            distance = pygplates.GeometryOnSphere.distance(
-                prev_point,
-                prev_boundary_geometry,
-                distance_threshold_radians=float(distance_threshold_radians),
-            )
-            if distance is not None:
-                # Detected a collision.
-                return True
-
-        return False
-
-
-_DEFAULT_COLLISION = _DefaultCollision()
-
-
-class _ContinentCollision(object):
-    """
-    Continental collision detection function class (the function is the '__call__' method).
-    """
-
-    def __init__(
-        self,
-        grd_output_dir,
-        chain_collision_detection=_DEFAULT_COLLISION,
-        verbose=False,
-    ):
-        """
-        grd_output_dir: The directory containing the continental grids.
-
-        chain_collision_detection: Another collision detection class/function to reference if we find no collision.
-                                   If None then no collision detection is chained. Defaults to the default collision detection.
-
-        verbose: Print progress messages
-        """
-
-        self.grd_output_dir = grd_output_dir
-        self.chain_collision_detection = chain_collision_detection
-        self.verbose = verbose
-
-        # Load a new grid each time the reconstruction time changes.
-        self.grid_time = None
-
-    @property
-    def grid_time(self):
-        return self._grid_time
-
-    @grid_time.setter
-    def grid_time(self, time):
-        from .grids import read_netcdf_grid
-
-        if time is None:
-            self._grid_time = time
-        else:
-            filename = "{:s}".format(self.grd_output_dir.format(time))
-            if self.verbose:
-                print(
-                    "Points masked against grid: {0}".format(os.path.basename(filename))
-                )
-            gridZ, gridX, gridY = read_netcdf_grid(filename, return_grids=True)
-            self.gridZ = gridZ
-            self.ni, self.nj = gridZ.shape
-            self.xmin = np.nanmin(gridX)
-            self.xmax = np.nanmax(gridX)
-            self.ymin = np.nanmin(gridY)
-            self.ymax = np.nanmax(gridY)
-            self._grid_time = float(time)
-
-    def __call__(
-        self,
-        rotation_model,
-        time,
-        reconstruction_time_interval,
-        prev_point,
-        curr_point,
-        prev_topology_plate_id,
-        prev_resolved_plate_boundary,
-        curr_topology_plate_id,
-        curr_resolved_plate_boundary,
-    ):
-        """
-        Returns True if a collision with a continent was detected, or returns result of
-        chained collision detection if 'self.chain_collision_detection' is not None.
-        """
-        # Load the grid for the current time if encountering a new time.
-        if time != self.grid_time:
-            self.grid_time = time
-            self.continent_deletion_count = 0
-
-        # Sample mask grid, which is one over continents and zero over oceans.
-        point_lat, point_lon = curr_point.to_lat_lon()
-        point_i = (self.ni - 1) * ((point_lat - self.ymin) / (self.ymax - self.ymin))
-        point_j = (self.nj - 1) * ((point_lon - self.xmin) / (self.xmax - self.xmin))
-        point_i_uint = np.rint(point_i).astype(np.uint)
-        point_j_uint = np.rint(point_j).astype(np.uint)
-        try:
-            mask_value = self.gridZ[point_i_uint, point_j_uint]
-        except IndexError:
-            point_i = np.clip(np.rint(point_i), 0, self.ni - 1).astype(np.int_)
-            point_j = np.clip(np.rint(point_j), 0, self.nj - 1).astype(np.int_)
-            mask_value = self.gridZ[point_i, point_j]
-        if mask_value >= 0.5:
-            # Detected a collision.
-            self.continent_deletion_count += 1
-            return True
-
-        # We didn't find a collision, so ask the chained collision detection if it did (if we have anything chained).
-        if self.chain_collision_detection:
-            return self.chain_collision_detection(
-                rotation_model,
-                time,
-                reconstruction_time_interval,
-                prev_point,
-                curr_point,
-                prev_topology_plate_id,
-                prev_resolved_plate_boundary,
-                curr_topology_plate_id,
-                curr_resolved_plate_boundary,
-            )
-
-        return False
-
-
-def reconstruct_points(
-    rotation_features_or_model,
-    topology_features,
-    reconstruction_begin_time,
-    reconstruction_end_time,
-    reconstruction_time_interval,
-    points,
-    point_begin_times=None,
-    point_end_times=None,
-    point_plate_ids=None,
-    detect_collisions=_DEFAULT_COLLISION,
-):
-    """
-    Function to reconstruct points using the ReconstructByTopologies class below.
-
-    For description of parameters see the ReconstructByTopologies class below.
-    """
-
-    topology_reconstruction = _ReconstructByTopologies(
-        rotation_features_or_model,
-        topology_features,
-        reconstruction_begin_time,
-        reconstruction_end_time,
-        reconstruction_time_interval,
-        points,
-        point_begin_times,
-        point_end_times,
-        point_plate_ids,
-        detect_collisions,
-    )
-
-    return topology_reconstruction.reconstruct()
-
-
-class _ReconstructByTopologies(object):
-    """
-    Class to reconstruct geometries using topologies.
-
-    Currently only points are supported.
-
-    use_plate_partitioner: If True then use pygplates.PlatePartitioner to partition points,
-                           otherwise use faster points_in_polygons.find_polygons().
-    """
-
-    use_plate_partitioner = False
-
-    def __init__(
-        self,
-        rotation_features_or_model,
-        topology_features,
-        reconstruction_begin_time,
-        reconstruction_end_time,
-        reconstruction_time_interval,
-        points,
-        point_begin_times=None,
-        point_end_times=None,
-        point_plate_ids=None,
-        detect_collisions=_DEFAULT_COLLISION,
-    ):
-        """
-        rotation_features_or_model: Rotation model or feature collection(s), or list of features, or filename(s).
-
-        topology_features: Topology feature collection(s), or list of features, or filename(s) or any combination of those.
-
-        detect_collisions: Collision detection function, or None. Defaults to _DEFAULT_COLLISION.
-        """
-
-        # Turn rotation data into a RotationModel (if not already).
-        if not isinstance(rotation_features_or_model, pygplates.RotationModel):
-            rotation_model = pygplates.RotationModel(rotation_features_or_model)
-        else:
-            rotation_model = rotation_features_or_model
-        self.rotation_model = rotation_model
-
-        # Turn topology data into a list of features (if not already).
-        self.topology_features = pygplates.FeaturesFunctionArgument(
-            topology_features
-        ).get_features()
-
-        # Set up an array of reconstruction times covering the reconstruction time span.
-        self.reconstruction_begin_time = reconstruction_begin_time
-        self.reconstruction_end_time = reconstruction_end_time
-        if reconstruction_time_interval <= 0.0:
-            raise ValueError("'reconstruction_time_interval' must be positive.")
-        # Reconstruction can go forward or backward in time.
-        if self.reconstruction_begin_time > self.reconstruction_end_time:
-            self.reconstruction_time_step = -reconstruction_time_interval
-        else:
-            self.reconstruction_time_step = reconstruction_time_interval
-        # Get number of times including end time if time span is a multiple of time step.
-        # The '1' is because, for example, 2 time intervals is 3 times.
-        # The '1e-6' deals with limited floating-point precision, eg, we want (3.0 - 0.0) / 1.0 to be 3.0 and not 2.999999 (which gets truncated to 2).
-        self.num_times = 1 + int(
-            math.floor(
-                1e-6
-                + float(self.reconstruction_end_time - self.reconstruction_begin_time)
-                / self.reconstruction_time_step
-            )
-        )
-        # It's possible the time step is larger than the time span, in which case we change it to equal the time span.
-        # This guarantees there'll be at least one time step (which has two times; one at either end of interval).
-        if self.num_times == 1:
-            self.num_times = 2
-            self.reconstruction_time_step = (
-                self.reconstruction_end_time - self.reconstruction_begin_time
-            )
-        self.reconstruction_time_interval = math.fabs(self.reconstruction_time_step)
-
-        self.last_time_index = self.num_times - 1
-
-        self.points = points
-        self.num_points = len(points)
-
-        # Use the specified point begin times if provided (otherwise use 'inf').
-        self.point_begin_times = point_begin_times
-        if self.point_begin_times is None:
-            self.point_begin_times = [float("inf")] * self.num_points
-        elif len(self.point_begin_times) != self.num_points:
-            raise ValueError(
-                "Length of 'point_begin_times' must match length of 'points'."
-            )
-
-        # Use the specified point end times if provided (otherwise use '-inf').
-        self.point_end_times = point_end_times
-        if self.point_end_times is None:
-            self.point_end_times = [float("-inf")] * self.num_points
-        elif len(self.point_end_times) != self.num_points:
-            raise ValueError(
-                "Length of 'point_end_times' must match length of 'points'."
-            )
-
-        # Use the specified point plate IDs if provided (otherwise use '0').
-        # These plate IDs are only used when a point falls outside all resolved topologies during a time step.
-        self.point_plate_ids = point_plate_ids
-        if self.point_plate_ids is None:
-            self.point_plate_ids = [0] * self.num_points
-        elif len(self.point_plate_ids) != self.num_points:
-            raise ValueError(
-                "Length of 'point_plate_ids' must match length of 'points'."
-            )
-
-        self.detect_collisions = detect_collisions
-
-    def reconstruct(self):
-        # Initialise the reconstruction.
-        self.begin_reconstruction()
-
-        # Loop over the reconstruction times until reached end of the reconstruction time span, or
-        # all points have entered their valid time range *and* either exited their time range or
-        # have been deactivated (subducted forward in time or consumed by MOR backward in time).
-        while self.reconstruct_to_next_time():
-            pass
-
-        return self.get_active_current_points()
-
-    def begin_reconstruction(self):
-        self.current_time_index = 0
-
-        # Set up point arrays.
-        # Store active and inactive points here (inactive points have None in corresponding entries).
-        self.prev_points = [None] * self.num_points
-        self.curr_points = [None] * self.num_points
-        self.next_points = [None] * self.num_points
-
-        # Each point can only get activated once (after deactivation it cannot be reactivated).
-        self.point_has_been_activated = [False] * self.num_points
-        self.num_activated_points = 0
-
-        # Set up topology arrays (corresponding to active/inactive points at same indices).
-        self.prev_topology_plate_ids = [None] * self.num_points
-        self.curr_topology_plate_ids = [None] * self.num_points
-        self.prev_resolved_plate_boundaries = [None] * self.num_points
-        self.curr_resolved_plate_boundaries = [None] * self.num_points
-
-        # Array to store indices of points found in continents
-        self.in_continent_indices = [None] * self.num_points
-        self.in_continent_points = [None] * self.num_points
-
-        self.deletedpoints = []
-
-        self._activate_deactivate_points()
-        self._find_resolved_topologies_containing_points()
-
-    def get_current_time(self):
-        return (
-            self.reconstruction_begin_time
-            + self.current_time_index * self.reconstruction_time_step
-        )
-
-    def get_all_current_points(self):
-        return self.curr_points
-
-    def get_active_current_points(self):
-        # Return only the active points (the ones that are not None).
-        return [point for point in self.get_all_current_points() if point is not None]
-
-    def get_in_continent_indices(self):
-        return self.in_continent_points, self.in_continent_indices
-
-    def reconstruct_to_next_time(self):
-        # If we're at the last time then there is no next time to reconstruct to.
-        if self.current_time_index == self.last_time_index:
-            return False
-
-        # If all points have been previously activated, but none are currently active then we're finished.
-        # This means all points have entered their valid time range *and* either exited their time range or
-        # have been deactivated (subducted forward in time or consumed by MOR backward in time).
-        if self.num_activated_points == self.num_points and not any(self.curr_points):
-            return False
-
-        # Cache stage rotations by plate ID.
-        # Use different dicts since using different rotation models and time steps, etc.
-        reconstruct_stage_rotation_dict = {}
-
-        current_time = self.get_current_time()
-
-        # Iterate over all points to reconstruct them to the next time step.
-        for point_index in range(self.num_points):
-            curr_point = self.curr_points[point_index]
-            if curr_point is None:
-                # Current point is not currently active.
-                # So we cannot reconstruct to next time.
-                self.next_points[point_index] = None
-                continue
-
-            # Get plate ID of resolved topology containing current point
-            # (this was determined in last call to '_find_resolved_topologies_containing_points()').
-            curr_plate_id = self.curr_topology_plate_ids[point_index]
-            if curr_plate_id is None:
-                # Current point is currently active but it fell outside all resolved polygons.
-                # So instead we just reconstruct using its plate ID (that was manually assigned by the user/caller).
-                curr_plate_id = self.point_plate_ids[point_index]
-
-            # Get the stage rotation that will move the point from where it is at the current time to its
-            # location at the next time step, based on the plate id that contains the point at the current time.
-
-            # Speed up by caching stage rotations in a dict.
-            stage_rotation = reconstruct_stage_rotation_dict.get(curr_plate_id)
-            if not stage_rotation:
-                stage_rotation = self.rotation_model.get_rotation(
-                    # Positive/negative time step means reconstructing backward/forward in time.
-                    current_time + self.reconstruction_time_step,
-                    curr_plate_id,
-                    current_time,
-                )
-                reconstruct_stage_rotation_dict[curr_plate_id] = stage_rotation
-
-            # Use the stage rotation to reconstruct the tracked point from position at current time
-            # to position at the next time step.
-            self.next_points[point_index] = stage_rotation * curr_point
-
-        #
-        # Set up for next loop iteration.
-        #
-        # Rotate previous, current and next point arrays.
-        # The new previous will be the old current.
-        # The new current will be the old next.
-        # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory).
-        self.prev_points, self.curr_points, self.next_points = (
-            self.curr_points,
-            self.next_points,
-            self.prev_points,
-        )
-        # Swap previous and current topology arrays.
-        # The new previous will be the old current.
-        # The new current will be the old previous (but values are ignored and overridden in next time step; just re-using its memory).
-        self.prev_topology_plate_ids, self.curr_topology_plate_ids = (
-            self.curr_topology_plate_ids,
-            self.prev_topology_plate_ids,
-        )
-        self.prev_resolved_plate_boundaries, self.curr_resolved_plate_boundaries = (
-            self.curr_resolved_plate_boundaries,
-            self.prev_resolved_plate_boundaries,
-        )
-
-        # Move the current time to the next time.
-        self.current_time_index += 1
-        current_time = self.get_current_time()
-
-        self._activate_deactivate_points()
-        self._find_resolved_topologies_containing_points()
-
-        # Iterate over all points to detect collisions.
-        if self.detect_collisions:
-            for point_index in range(self.num_points):
-                prev_point = self.prev_points[point_index]
-                curr_point = self.curr_points[point_index]
-                if prev_point is None or curr_point is None:
-                    # If current point is not currently active then no need to detect a collision for it (to deactivate it).
-                    # Also previous point might just have been activated now, at end of current time step, and hence
-                    # not active at beginning of time step.
-                    continue
-
-                # Get plate IDs of resolved topology containing previous and current point
-                # (this was determined in last call to '_find_resolved_topologies_containing_points()').
-                #
-                # Note that could be None, so the collision detection needs to handle that.
-                prev_plate_id = self.prev_topology_plate_ids[point_index]
-                curr_plate_id = self.curr_topology_plate_ids[point_index]
-
-                # Detect collisions at the end of the current time step since we need previous, and current, points and topologies.
-                # De-activate point (in 'curr_points') if subducted (forward in time) or consumed back into MOR (backward in time).
-                if self.detect_collisions(
-                    self.rotation_model,
-                    current_time,
-                    self.reconstruction_time_interval,
-                    prev_point,
-                    curr_point,
-                    prev_plate_id,
-                    self.prev_resolved_plate_boundaries[point_index],
-                    curr_plate_id,
-                    self.curr_resolved_plate_boundaries[point_index],
-                ):
-                    # An inactive point in 'curr_points' becomes None.
-                    # It may have been reconstructed from the previous time step to a valid position
-                    # but now we override that result as inactive.
-                    self.curr_points[point_index] = None
-                    # self.curr_points.remove(self.curr_points[point_index])
-                    self.deletedpoints.append(point_index)
-
-        # We successfully reconstructed to the next time.
-        return True
-
-    def _activate_deactivate_points(self):
-        current_time = self.get_current_time()
-
-        # Iterate over all points and activate/deactivate as necessary depending on each point's valid time range.
-        for point_index in range(self.num_points):
-            if self.curr_points[point_index] is None:
-                if not self.point_has_been_activated[point_index]:
-                    # Point is not active and has never been activated, so see if can activate it.
-                    if (
-                        current_time <= self.point_begin_times[point_index]
-                        and current_time >= self.point_end_times[point_index]
-                    ):
-                        # The initial point is assumed to be the position at the current time
-                        # which is typically the point's begin time (approximately).
-                        # But it could be the beginning of the reconstruction time span (specified in constructor)
-                        # if that falls in the middle of the point's valid time range - in this case the
-                        # initial point position is assumed to be in a position that is some time *after*
-                        # it appeared (at its begin time) - and this can happen, for example, if you have a
-                        # uniform grids of points at some intermediate time and want to see how they
-                        # reconstruct to either a younger or older time (remembering that points can
-                        # be subducted forward in time and consumed back into a mid-ocean ridge going
-                        # backward in time).
-                        self.curr_points[point_index] = self.points[point_index]
-                        self.point_has_been_activated[point_index] = True
-                        self.num_activated_points += 1
-            else:
-                # Point is active, so see if can deactivate it.
-                if not (
-                    current_time <= self.point_begin_times[point_index]
-                    and current_time >= self.point_end_times[point_index]
-                ):
-                    self.curr_points[point_index] = None
-
-    def _find_resolved_topologies_containing_points(self):
-        from . import ptt as _ptt
-
-        current_time = self.get_current_time()
-
-        # Resolve the plate polygons for the current time.
-        resolved_topologies = []
-        pygplates.resolve_topologies(
-            self.topology_features,
-            self.rotation_model,
-            resolved_topologies,
-            current_time,
-        )
-
-        if _ReconstructByTopologies.use_plate_partitioner:
-            # Create a plate partitioner from the resolved polygons.
-            plate_partitioner = pygplates.PlatePartitioner(
-                resolved_topologies, self.rotation_model
-            )
-        else:
-            # Some of 'curr_points' will be None so 'curr_valid_points' contains only the valid (not None)
-            # points, and 'curr_valid_points_indices' is the same length as 'curr_points' but indexes into
-            # 'curr_valid_points' so we can quickly find which point (and hence which resolved topology)
-            # in 'curr_valid_points' is associated with the a particular point in 'curr_points'.
-            curr_valid_points = []
-            curr_valid_points_indices = [None] * self.num_points
-            for point_index, curr_point in enumerate(self.curr_points):
-                if curr_point is not None:
-                    curr_valid_points_indices[point_index] = len(curr_valid_points)
-                    curr_valid_points.append(curr_point)
-            # For each valid current point find the resolved topology containing it.
-            resolved_topologies_containing_curr_valid_points = (
-                _ptt.utils.points_in_polygons.find_polygons(
-                    curr_valid_points,
-                    [
-                        resolved_topology.get_resolved_boundary()
-                        for resolved_topology in resolved_topologies
-                    ],
-                    resolved_topologies,
-                )
-            )
-
-        # Iterate over all points.
-        for point_index, curr_point in enumerate(self.curr_points):
-            if curr_point is None:
-                # Current point is not currently active - so skip it.
-                self.curr_topology_plate_ids[point_index] = None
-                self.curr_resolved_plate_boundaries[point_index] = None
-                continue
-
-            # Find the plate id of the polygon that contains 'curr_point'.
-            if _ReconstructByTopologies.use_plate_partitioner:
-                curr_polygon = plate_partitioner.partition_point(curr_point)
-            else:
-                curr_polygon = resolved_topologies_containing_curr_valid_points[
-                    # Index back into 'curr_valid_points' and hence also into
-                    # 'resolved_topologies_containing_curr_valid_points'.
-                    curr_valid_points_indices[point_index]
-                ]
-            self.curr_resolved_plate_boundaries[point_index] = curr_polygon
-
-            # If the polygon is None, that means (presumably) that it fell into a crack between
-            # topologies. So it will be skipped and thrown away from future iterations.
-            if curr_polygon is None:
-                self.curr_topology_plate_ids[point_index] = None
-                continue
-
-            # Set the plate ID of resolved topology containing current point.
-            self.curr_topology_plate_ids[point_index] = (
-                curr_polygon.get_feature().get_reconstruction_plate_id()
-            )
diff --git a/gplately/utils/__init__.py b/gplately/utils/__init__.py
index 301f4199..6dd0a23a 100644
--- a/gplately/utils/__init__.py
+++ b/gplately/utils/__init__.py
@@ -22,6 +22,10 @@
     create_icosahedral_mesh,
     ensure_polygon_geometry,
     point_in_polygon_routine,
+    DefaultCollision,
+    ContinentCollision,
+    reconstruct_points,
+    ReconstructByTopologies,
 )
 
 __all__ = [
@@ -32,6 +36,10 @@
     "create_icosahedral_mesh",
     "ensure_polygon_geometry",
     "point_in_polygon_routine",
+    "DefaultCollision",
+    "ContinentCollision",
+    "econstruct_points",
+    "ReconstructByTopologies",
 ]
 
 __pdoc__ = {
diff --git a/gplately/utils/seafloor_grid_utils.py b/gplately/utils/seafloor_grid_utils.py
index 5d68e8f0..d196aff8 100644
--- a/gplately/utils/seafloor_grid_utils.py
+++ b/gplately/utils/seafloor_grid_utils.py
@@ -18,10 +18,13 @@
 Auxiliary functions for SeafloorGrid
 """
 
+import math
+import os
+
 import numpy as np
 import pygplates
 
-from .. import ptt
+from .. import ptt as _ptt
 
 
 def create_icosahedral_mesh(refinement_levels):
@@ -189,7 +192,7 @@ def point_in_polygon_routine(multi_point, COB_polygons):
         polygons[ind] = geom
     proxies = np.ones(polygons.size)
 
-    pip_result = ptt.utils.points_in_polygons.find_polygons(
+    pip_result = _ptt.utils.points_in_polygons.find_polygons(
         multi_point, polygons, proxies, all_polygons=False
     )  # 1 for points in polygons, None for points outside
     zvals = np.array(
@@ -206,3 +209,783 @@ def point_in_polygon_routine(multi_point, COB_polygons):
         pygplates.MultiPointOnSphere(points_out_arr),
         zvals,
     )
+
+# FROM RECONSTRUCT_BY_TOPOLOGIES.PY
+class DefaultCollision(object):
+    """
+    Default collision detection function class (the function is the '__call__' method).
+    """
+
+    DEFAULT_GLOBAL_COLLISION_PARAMETERS = (7.0, 10.0)
+    """
+    Default collision parameters for all feature types.
+
+    This is a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my):
+    Here we default to the same constants used internally in GPlates 2.0 (ie, 7.0 and 10.0).
+    """
+
+    def __init__(
+        self,
+        global_collision_parameters=DEFAULT_GLOBAL_COLLISION_PARAMETERS,
+        feature_specific_collision_parameters=None,
+    ):
+        """
+        global_collision_parameters: The collision parameters to use for any feature type not specified in 'feature_specific_collision_parameters'.
+                                     Should be a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my).
+                                     The first threshold parameter means:
+                                        A point that transitions from one plate to another can disappear if the change in velocity exceeds this threshold.
+                                     The second threshold parameter means:
+                                        Only those transitioning points exceeding the threshold velocity delta and that are close enough to a plate boundary can disappear.
+                                        The distance is proportional to the relative velocity (change in velocity), plus a constant offset based on the threshold distance to boundary
+                                        to account for plate boundaries that change shape significantly from one time step to the next
+                                        (note that some boundaries are meant to do this and others are a result of digitisation).
+                                        The actual distance threshold used is (threshold_distance_to_boundary + relative_velocity) * time_interval
+                                     Defaults to parameters used in GPlates 2.0, if not specified.
+
+        feature_specific_collision_parameters: Optional sequence of collision parameters specific to feature types.
+                                               If specified then should be a sequence of 2-tuples, with each 2-tuple specifying (feature_type, collision_parameters).
+                                               And where each 'collision_parameters' is a 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my).
+                                                   See 'global_collision_parameters' for details on these thresholds.
+                                               Any feature type not specified here defaults to using 'global_collision_parameters'.
+        """
+
+        # Convert list of (feature_type, collision_parameters) tuples to a dictionary.
+        if feature_specific_collision_parameters:
+            self.feature_specific_collision_parameters = dict(
+                feature_specific_collision_parameters
+            )
+        else:
+            self.feature_specific_collision_parameters = dict()
+        # Fallback for any feature type not specified in the optional feature-specific list.
+        self.global_collision_parameters = global_collision_parameters
+
+        # Used to improve performance by caching velocity stage rotations in a dict (for a specific reconstruction time).
+        self.velocity_stage_rotation_dict = {}
+        self.velocity_stage_rotation_time = None
+
+    def __call__(
+        self,
+        rotation_model,
+        time,
+        reconstruction_time_interval,
+        prev_point,
+        curr_point,
+        prev_topology_plate_id,
+        prev_resolved_plate_boundary,
+        curr_topology_plate_id,
+        curr_resolved_plate_boundary,
+    ):
+        """
+        Returns True if a collision was detected.
+
+        If transitioning from a rigid plate to another rigid plate with a different plate ID then
+        calculate the difference in velocities and continue testing as follows
+        (otherwise, if there's no transition, then the point is still active)...
+
+        If the velocity difference is below a threshold then we assume the previous plate was split,
+        or two plates joined. In this case the point has not subducted (forward in time) or
+        been consumed by a mid-ocean (backward in time) and hence is still active.
+
+        If the velocity difference is large enough then we see if the distance of the *previous* position
+        to the polygon boundary (of rigid plate containing it) exceeds a threshold.
+        If the distance exceeds the threshold then the point is far enough away from the boundary that it
+        cannot be subducted or consumed by it and hence the point is still active.
+        However if the point is close enough then we assume the point was subducted/consumed
+        (remember that the point switched plate IDs).
+        Also note that the threshold distance increases according to the velocity difference to account for fast
+        moving points (that would otherwise tunnel through the boundary and accrete onto the other plate).
+        The reason for testing the distance from the *previous* point, and not from the *current* point, is:
+
+          (i)  A topological boundary may *appear* near the current point (such as a plate split at the current time)
+               and we don't want that split to consume the current point regardless of the velocity difference.
+               It won't get consumed because the *previous* point was not near a boundary (because before split happened).
+               If the velocity difference is large enough then it might cause the current point to transition to the
+               adjacent split plate in the *next* time step (and that's when it should get consumed, not in the current time step).
+               An example of this is a mid-ocean ridge suddenly appearing (going forward in time).
+
+          (ii) A topological boundary may *disappear* near the current point (such as a plate merge at the current time)
+               and we want that merge to consume the current point if the velocity difference is large enough.
+               In this case the *previous* point is near a boundary (because before plate merged) and hence can be
+               consumed (provided velocity difference is large enough). And since the boundary existed in the previous
+               time step, it will affect position of the current point (and whether it gets consumed or not).
+               An example of this is a mid-ocean ridge suddenly disappearing (going backward in time).
+
+        ...note that items (i) and (ii) above apply both going forward and backward in time.
+        """
+
+        # See if a collision occurred.
+        if (
+            curr_topology_plate_id != prev_topology_plate_id
+            and prev_topology_plate_id is not None
+            and curr_topology_plate_id is not None
+        ):
+            #
+            # Speed up by caching velocity stage rotations in a dict.
+            #
+            if time != self.velocity_stage_rotation_time:
+                # We've just switched to a new time so clear the cache.
+                #
+                # We only cache stage rotations for a specific time.
+                # We only really need to cache different plate IDs at the same 'time', so this avoids caching for all times
+                # (which would also require including 'time' in the key) and using memory unnecessarily.
+                self.velocity_stage_rotation_dict.clear()
+                self.velocity_stage_rotation_time = time
+            prev_location_velocity_stage_rotation = (
+                self.velocity_stage_rotation_dict.get(prev_topology_plate_id)
+            )
+            if not prev_location_velocity_stage_rotation:
+                prev_location_velocity_stage_rotation = rotation_model.get_rotation(
+                    time + 1, prev_topology_plate_id, time
+                )
+                self.velocity_stage_rotation_dict[prev_topology_plate_id] = (
+                    prev_location_velocity_stage_rotation
+                )
+            curr_location_velocity_stage_rotation = (
+                self.velocity_stage_rotation_dict.get(curr_topology_plate_id)
+            )
+            if not curr_location_velocity_stage_rotation:
+                curr_location_velocity_stage_rotation = rotation_model.get_rotation(
+                    time + 1, curr_topology_plate_id, time
+                )
+                self.velocity_stage_rotation_dict[curr_topology_plate_id] = (
+                    curr_location_velocity_stage_rotation
+                )
+
+            # Note that even though the current point is not inside the previous boundary (because different plate ID), we can still
+            # calculate a velocity using its plate ID (because we really should use the same point in our velocity comparison).
+            prev_location_velocity = pygplates.calculate_velocities(
+                (curr_point,),
+                prev_location_velocity_stage_rotation,
+                1,
+                pygplates.VelocityUnits.kms_per_my,
+            )[0]
+            curr_location_velocity = pygplates.calculate_velocities(
+                (curr_point,),
+                curr_location_velocity_stage_rotation,
+                1,
+                pygplates.VelocityUnits.kms_per_my,
+            )[0]
+
+            delta_velocity = curr_location_velocity - prev_location_velocity
+            delta_velocity_magnitude = delta_velocity.get_magnitude()
+
+            # If we have feature-specific collision parameters then iterate over the boundary sub-segments of the *previous* topological boundary
+            # and test proximity to each sub-segment individually (with sub-segment feature type specific collision parameters).
+            # Otherwise just test proximity to the entire boundary polygon using the global collision parameters.
+            if self.feature_specific_collision_parameters:
+                for (
+                    prev_boundary_sub_segment
+                ) in prev_resolved_plate_boundary.get_boundary_sub_segments():
+                    # Use feature-specific collision parameters if found (falling back to global collision parameters).
+                    (
+                        threshold_velocity_delta,
+                        threshold_distance_to_boundary_per_my,
+                    ) = self.feature_specific_collision_parameters.get(
+                        prev_boundary_sub_segment.get_feature().get_feature_type(),
+                        # Default to global collision parameters if no collision parameters specified for sub-segment's feature type...
+                        self.global_collision_parameters,
+                    )
+
+                    # Since each feature type could use different collision parameters we must use the current boundary sub-segment instead of the boundary polygon.
+                    if self._detect_collision_using_collision_parameters(
+                        reconstruction_time_interval,
+                        delta_velocity_magnitude,
+                        prev_point,
+                        prev_boundary_sub_segment.get_resolved_geometry(),
+                        threshold_velocity_delta,
+                        threshold_distance_to_boundary_per_my,
+                    ):
+                        # Detected a collision.
+                        return True
+            else:
+                # No feature-specific collision parameters so use global fallback.
+                (
+                    threshold_velocity_delta,
+                    threshold_distance_to_boundary_per_my,
+                ) = self.global_collision_parameters
+
+                # Since all feature types use the same collision parameters we can use the boundary polygon instead of iterating over its sub-segments.
+                if self._detect_collision_using_collision_parameters(
+                    reconstruction_time_interval,
+                    delta_velocity_magnitude,
+                    prev_point,
+                    prev_resolved_plate_boundary.get_resolved_boundary(),
+                    threshold_velocity_delta,
+                    threshold_distance_to_boundary_per_my,
+                ):
+                    # Detected a collision.
+                    return True
+
+        return False
+
+    def _detect_collision_using_collision_parameters(
+        self,
+        reconstruction_time_interval,
+        delta_velocity_magnitude,
+        prev_point,
+        prev_boundary_geometry,
+        threshold_velocity_delta,
+        threshold_distance_to_boundary_per_my,
+    ):
+        if delta_velocity_magnitude > threshold_velocity_delta:
+            # Add the minimum distance threshold to the delta velocity threshold.
+            # The delta velocity threshold only allows those points that are close enough to the boundary to reach
+            # it given their current relative velocity.
+            # The minimum distance threshold accounts for sudden changes in the shape of a plate boundary
+            # which are no supposed to represent a new or shifted boundary but are just a result of the topology
+            # builder/user digitising a new boundary line that differs noticeably from that of the previous time period.
+            distance_threshold_radians = (
+                (threshold_distance_to_boundary_per_my + delta_velocity_magnitude)
+                * reconstruction_time_interval
+                / pygplates.Earth.equatorial_radius_in_kms
+            )
+            distance_threshold_radians = min(distance_threshold_radians, math.pi)
+            distance_threshold_radians = max(distance_threshold_radians, 0.0)
+
+            distance = pygplates.GeometryOnSphere.distance(
+                prev_point,
+                prev_boundary_geometry,
+                distance_threshold_radians=float(distance_threshold_radians),
+            )
+            if distance is not None:
+                # Detected a collision.
+                return True
+
+        return False
+
+
+DEFAULT_COLLISION = DefaultCollision()
+
+
+class ContinentCollision(object):
+    """
+    Continental collision detection function class (the function is the '__call__' method).
+    """
+
+    def __init__(
+        self,
+        grd_output_dir,
+        chain_collision_detection=DEFAULT_COLLISION,
+        verbose=False,
+    ):
+        """
+        grd_output_dir: The directory containing the continental grids.
+
+        chain_collision_detection: Another collision detection class/function to reference if we find no collision.
+                                   If None then no collision detection is chained. Defaults to the default collision detection.
+
+        verbose: Print progress messages
+        """
+
+        self.grd_output_dir = grd_output_dir
+        self.chain_collision_detection = chain_collision_detection
+        self.verbose = verbose
+
+        # Load a new grid each time the reconstruction time changes.
+        self.grid_time = None
+
+    @property
+    def grid_time(self):
+        return self._grid_time
+
+    @grid_time.setter
+    def grid_time(self, time):
+        from ..grids import read_netcdf_grid
+
+        if time is None:
+            self._grid_time = time
+        else:
+            filename = "{:s}".format(self.grd_output_dir.format(time))
+            if self.verbose:
+                print(
+                    "Points masked against grid: {0}".format(os.path.basename(filename))
+                )
+            gridZ, gridX, gridY = read_netcdf_grid(filename, return_grids=True)
+            self.gridZ = gridZ
+            self.ni, self.nj = gridZ.shape
+            self.xmin = np.nanmin(gridX)
+            self.xmax = np.nanmax(gridX)
+            self.ymin = np.nanmin(gridY)
+            self.ymax = np.nanmax(gridY)
+            self._grid_time = float(time)
+
+    def __call__(
+        self,
+        rotation_model,
+        time,
+        reconstruction_time_interval,
+        prev_point,
+        curr_point,
+        prev_topology_plate_id,
+        prev_resolved_plate_boundary,
+        curr_topology_plate_id,
+        curr_resolved_plate_boundary,
+    ):
+        """
+        Returns True if a collision with a continent was detected, or returns result of
+        chained collision detection if 'self.chain_collision_detection' is not None.
+        """
+        # Load the grid for the current time if encountering a new time.
+        if time != self.grid_time:
+            self.grid_time = time
+            self.continent_deletion_count = 0
+
+        # Sample mask grid, which is one over continents and zero over oceans.
+        point_lat, point_lon = curr_point.to_lat_lon()
+        point_i = (self.ni - 1) * ((point_lat - self.ymin) / (self.ymax - self.ymin))
+        point_j = (self.nj - 1) * ((point_lon - self.xmin) / (self.xmax - self.xmin))
+        point_i_uint = np.rint(point_i).astype(np.uint)
+        point_j_uint = np.rint(point_j).astype(np.uint)
+        try:
+            mask_value = self.gridZ[point_i_uint, point_j_uint]
+        except IndexError:
+            point_i = np.clip(np.rint(point_i), 0, self.ni - 1).astype(np.int_)
+            point_j = np.clip(np.rint(point_j), 0, self.nj - 1).astype(np.int_)
+            mask_value = self.gridZ[point_i, point_j]
+        if mask_value >= 0.5:
+            # Detected a collision.
+            self.continent_deletion_count += 1
+            return True
+
+        # We didn't find a collision, so ask the chained collision detection if it did (if we have anything chained).
+        if self.chain_collision_detection:
+            return self.chain_collision_detection(
+                rotation_model,
+                time,
+                reconstruction_time_interval,
+                prev_point,
+                curr_point,
+                prev_topology_plate_id,
+                prev_resolved_plate_boundary,
+                curr_topology_plate_id,
+                curr_resolved_plate_boundary,
+            )
+
+        return False
+
+
+def reconstruct_points(
+    rotation_features_or_model,
+    topology_features,
+    reconstruction_begin_time,
+    reconstruction_end_time,
+    reconstruction_time_interval,
+    points,
+    point_begin_times=None,
+    point_end_times=None,
+    point_plate_ids=None,
+    detect_collisions=DEFAULT_COLLISION,
+):
+    """
+    Function to reconstruct points using the ReconstructByTopologies class below.
+
+    For description of parameters see the ReconstructByTopologies class below.
+    """
+
+    topology_reconstruction = ReconstructByTopologies(
+        rotation_features_or_model,
+        topology_features,
+        reconstruction_begin_time,
+        reconstruction_end_time,
+        reconstruction_time_interval,
+        points,
+        point_begin_times,
+        point_end_times,
+        point_plate_ids,
+        detect_collisions,
+    )
+
+    return topology_reconstruction.reconstruct()
+
+class ReconstructByTopologies(object):
+    """
+    Class to reconstruct geometries using topologies.
+
+    Currently only points are supported.
+
+    use_plate_partitioner: If True then use pygplates.PlatePartitioner to partition points,
+                           otherwise use faster points_in_polygons.find_polygons().
+    """
+
+    use_plate_partitioner = False
+
+    def __init__(
+        self,
+        rotation_features_or_model,
+        topology_features,
+        reconstruction_begin_time,
+        reconstruction_end_time,
+        reconstruction_time_interval,
+        points,
+        point_begin_times=None,
+        point_end_times=None,
+        point_plate_ids=None,
+        detect_collisions=DEFAULT_COLLISION,
+    ):
+        """
+        rotation_features_or_model: Rotation model or feature collection(s), or list of features, or filename(s).
+
+        topology_features: Topology feature collection(s), or list of features, or filename(s) or any combination of those.
+
+        detect_collisions: Collision detection function, or None. Defaults to DEFAULT_COLLISION.
+        """
+
+        # Turn rotation data into a RotationModel (if not already).
+        if not isinstance(rotation_features_or_model, pygplates.RotationModel):
+            rotation_model = pygplates.RotationModel(rotation_features_or_model)
+        else:
+            rotation_model = rotation_features_or_model
+        self.rotation_model = rotation_model
+
+        # Turn topology data into a list of features (if not already).
+        self.topology_features = pygplates.FeaturesFunctionArgument(
+            topology_features
+        ).get_features()
+
+        # Set up an array of reconstruction times covering the reconstruction time span.
+        self.reconstruction_begin_time = reconstruction_begin_time
+        self.reconstruction_end_time = reconstruction_end_time
+        if reconstruction_time_interval <= 0.0:
+            raise ValueError("'reconstruction_time_interval' must be positive.")
+        # Reconstruction can go forward or backward in time.
+        if self.reconstruction_begin_time > self.reconstruction_end_time:
+            self.reconstruction_time_step = -reconstruction_time_interval
+        else:
+            self.reconstruction_time_step = reconstruction_time_interval
+        # Get number of times including end time if time span is a multiple of time step.
+        # The '1' is because, for example, 2 time intervals is 3 times.
+        # The '1e-6' deals with limited floating-point precision, eg, we want (3.0 - 0.0) / 1.0 to be 3.0 and not 2.999999 (which gets truncated to 2).
+        self.num_times = 1 + int(
+            math.floor(
+                1e-6
+                + float(self.reconstruction_end_time - self.reconstruction_begin_time)
+                / self.reconstruction_time_step
+            )
+        )
+        # It's possible the time step is larger than the time span, in which case we change it to equal the time span.
+        # This guarantees there'll be at least one time step (which has two times; one at either end of interval).
+        if self.num_times == 1:
+            self.num_times = 2
+            self.reconstruction_time_step = (
+                self.reconstruction_end_time - self.reconstruction_begin_time
+            )
+        self.reconstruction_time_interval = math.fabs(self.reconstruction_time_step)
+
+        self.last_time_index = self.num_times - 1
+
+        self.points = points
+        self.num_points = len(points)
+
+        # Use the specified point begin times if provided (otherwise use 'inf').
+        self.point_begin_times = point_begin_times
+        if self.point_begin_times is None:
+            self.point_begin_times = [float("inf")] * self.num_points
+        elif len(self.point_begin_times) != self.num_points:
+            raise ValueError(
+                "Length of 'point_begin_times' must match length of 'points'."
+            )
+
+        # Use the specified point end times if provided (otherwise use '-inf').
+        self.point_end_times = point_end_times
+        if self.point_end_times is None:
+            self.point_end_times = [float("-inf")] * self.num_points
+        elif len(self.point_end_times) != self.num_points:
+            raise ValueError(
+                "Length of 'point_end_times' must match length of 'points'."
+            )
+
+        # Use the specified point plate IDs if provided (otherwise use '0').
+        # These plate IDs are only used when a point falls outside all resolved topologies during a time step.
+        self.point_plate_ids = point_plate_ids
+        if self.point_plate_ids is None:
+            self.point_plate_ids = [0] * self.num_points
+        elif len(self.point_plate_ids) != self.num_points:
+            raise ValueError(
+                "Length of 'point_plate_ids' must match length of 'points'."
+            )
+
+        self.detect_collisions = detect_collisions
+
+    def reconstruct(self):
+        # Initialise the reconstruction.
+        self.begin_reconstruction()
+
+        # Loop over the reconstruction times until reached end of the reconstruction time span, or
+        # all points have entered their valid time range *and* either exited their time range or
+        # have been deactivated (subducted forward in time or consumed by MOR backward in time).
+        while self.reconstruct_to_next_time():
+            pass
+
+        return self.get_active_current_points()
+
+    def begin_reconstruction(self):
+        self.current_time_index = 0
+
+        # Set up point arrays.
+        # Store active and inactive points here (inactive points have None in corresponding entries).
+        self.prev_points = [None] * self.num_points
+        self.curr_points = [None] * self.num_points
+        self.next_points = [None] * self.num_points
+
+        # Each point can only get activated once (after deactivation it cannot be reactivated).
+        self.point_has_been_activated = [False] * self.num_points
+        self.num_activated_points = 0
+
+        # Set up topology arrays (corresponding to active/inactive points at same indices).
+        self.prev_topology_plate_ids = [None] * self.num_points
+        self.curr_topology_plate_ids = [None] * self.num_points
+        self.prev_resolved_plate_boundaries = [None] * self.num_points
+        self.curr_resolved_plate_boundaries = [None] * self.num_points
+
+        # Array to store indices of points found in continents
+        self.in_continent_indices = [None] * self.num_points
+        self.in_continent_points = [None] * self.num_points
+
+        self.deletedpoints = []
+
+        self._activate_deactivate_points()
+        self._find_resolved_topologies_containing_points()
+
+    def get_current_time(self):
+        return (
+            self.reconstruction_begin_time
+            + self.current_time_index * self.reconstruction_time_step
+        )
+
+    def get_all_current_points(self):
+        return self.curr_points
+
+    def get_active_current_points(self):
+        # Return only the active points (the ones that are not None).
+        return [point for point in self.get_all_current_points() if point is not None]
+
+    def get_in_continent_indices(self):
+        return self.in_continent_points, self.in_continent_indices
+
+    def reconstruct_to_next_time(self):
+        # If we're at the last time then there is no next time to reconstruct to.
+        if self.current_time_index == self.last_time_index:
+            return False
+
+        # If all points have been previously activated, but none are currently active then we're finished.
+        # This means all points have entered their valid time range *and* either exited their time range or
+        # have been deactivated (subducted forward in time or consumed by MOR backward in time).
+        if self.num_activated_points == self.num_points and not any(self.curr_points):
+            return False
+
+        # Cache stage rotations by plate ID.
+        # Use different dicts since using different rotation models and time steps, etc.
+        reconstruct_stage_rotation_dict = {}
+
+        current_time = self.get_current_time()
+
+        # Iterate over all points to reconstruct them to the next time step.
+        for point_index in range(self.num_points):
+            curr_point = self.curr_points[point_index]
+            if curr_point is None:
+                # Current point is not currently active.
+                # So we cannot reconstruct to next time.
+                self.next_points[point_index] = None
+                continue
+
+            # Get plate ID of resolved topology containing current point
+            # (this was determined in last call to '_find_resolved_topologies_containing_points()').
+            curr_plate_id = self.curr_topology_plate_ids[point_index]
+            if curr_plate_id is None:
+                # Current point is currently active but it fell outside all resolved polygons.
+                # So instead we just reconstruct using its plate ID (that was manually assigned by the user/caller).
+                curr_plate_id = self.point_plate_ids[point_index]
+
+            # Get the stage rotation that will move the point from where it is at the current time to its
+            # location at the next time step, based on the plate id that contains the point at the current time.
+
+            # Speed up by caching stage rotations in a dict.
+            stage_rotation = reconstruct_stage_rotation_dict.get(curr_plate_id)
+            if not stage_rotation:
+                stage_rotation = self.rotation_model.get_rotation(
+                    # Positive/negative time step means reconstructing backward/forward in time.
+                    current_time + self.reconstruction_time_step,
+                    curr_plate_id,
+                    current_time,
+                )
+                reconstruct_stage_rotation_dict[curr_plate_id] = stage_rotation
+
+            # Use the stage rotation to reconstruct the tracked point from position at current time
+            # to position at the next time step.
+            self.next_points[point_index] = stage_rotation * curr_point
+
+        #
+        # Set up for next loop iteration.
+        #
+        # Rotate previous, current and next point arrays.
+        # The new previous will be the old current.
+        # The new current will be the old next.
+        # The new next will be the old previous (but values are ignored and overridden in next time step; just re-using its memory).
+        self.prev_points, self.curr_points, self.next_points = (
+            self.curr_points,
+            self.next_points,
+            self.prev_points,
+        )
+        # Swap previous and current topology arrays.
+        # The new previous will be the old current.
+        # The new current will be the old previous (but values are ignored and overridden in next time step; just re-using its memory).
+        self.prev_topology_plate_ids, self.curr_topology_plate_ids = (
+            self.curr_topology_plate_ids,
+            self.prev_topology_plate_ids,
+        )
+        self.prev_resolved_plate_boundaries, self.curr_resolved_plate_boundaries = (
+            self.curr_resolved_plate_boundaries,
+            self.prev_resolved_plate_boundaries,
+        )
+
+        # Move the current time to the next time.
+        self.current_time_index += 1
+        current_time = self.get_current_time()
+
+        self._activate_deactivate_points()
+        self._find_resolved_topologies_containing_points()
+
+        # Iterate over all points to detect collisions.
+        if self.detect_collisions:
+            for point_index in range(self.num_points):
+                prev_point = self.prev_points[point_index]
+                curr_point = self.curr_points[point_index]
+                if prev_point is None or curr_point is None:
+                    # If current point is not currently active then no need to detect a collision for it (to deactivate it).
+                    # Also previous point might just have been activated now, at end of current time step, and hence
+                    # not active at beginning of time step.
+                    continue
+
+                # Get plate IDs of resolved topology containing previous and current point
+                # (this was determined in last call to '_find_resolved_topologies_containing_points()').
+                #
+                # Note that could be None, so the collision detection needs to handle that.
+                prev_plate_id = self.prev_topology_plate_ids[point_index]
+                curr_plate_id = self.curr_topology_plate_ids[point_index]
+
+                # Detect collisions at the end of the current time step since we need previous, and current, points and topologies.
+                # De-activate point (in 'curr_points') if subducted (forward in time) or consumed back into MOR (backward in time).
+                if self.detect_collisions(
+                    self.rotation_model,
+                    current_time,
+                    self.reconstruction_time_interval,
+                    prev_point,
+                    curr_point,
+                    prev_plate_id,
+                    self.prev_resolved_plate_boundaries[point_index],
+                    curr_plate_id,
+                    self.curr_resolved_plate_boundaries[point_index],
+                ):
+                    # An inactive point in 'curr_points' becomes None.
+                    # It may have been reconstructed from the previous time step to a valid position
+                    # but now we override that result as inactive.
+                    self.curr_points[point_index] = None
+                    # self.curr_points.remove(self.curr_points[point_index])
+                    self.deletedpoints.append(point_index)
+
+        # We successfully reconstructed to the next time.
+        return True
+
+    def _activate_deactivate_points(self):
+        current_time = self.get_current_time()
+
+        # Iterate over all points and activate/deactivate as necessary depending on each point's valid time range.
+        for point_index in range(self.num_points):
+            if self.curr_points[point_index] is None:
+                if not self.point_has_been_activated[point_index]:
+                    # Point is not active and has never been activated, so see if can activate it.
+                    if (
+                        current_time <= self.point_begin_times[point_index]
+                        and current_time >= self.point_end_times[point_index]
+                    ):
+                        # The initial point is assumed to be the position at the current time
+                        # which is typically the point's begin time (approximately).
+                        # But it could be the beginning of the reconstruction time span (specified in constructor)
+                        # if that falls in the middle of the point's valid time range - in this case the
+                        # initial point position is assumed to be in a position that is some time *after*
+                        # it appeared (at its begin time) - and this can happen, for example, if you have a
+                        # uniform grids of points at some intermediate time and want to see how they
+                        # reconstruct to either a younger or older time (remembering that points can
+                        # be subducted forward in time and consumed back into a mid-ocean ridge going
+                        # backward in time).
+                        self.curr_points[point_index] = self.points[point_index]
+                        self.point_has_been_activated[point_index] = True
+                        self.num_activated_points += 1
+            else:
+                # Point is active, so see if can deactivate it.
+                if not (
+                    current_time <= self.point_begin_times[point_index]
+                    and current_time >= self.point_end_times[point_index]
+                ):
+                    self.curr_points[point_index] = None
+
+    def _find_resolved_topologies_containing_points(self):
+
+        current_time = self.get_current_time()
+
+        # Resolve the plate polygons for the current time.
+        resolved_topologies = []
+        pygplates.resolve_topologies(
+            self.topology_features,
+            self.rotation_model,
+            resolved_topologies,
+            current_time,
+        )
+
+        if ReconstructByTopologies.use_plate_partitioner:
+            # Create a plate partitioner from the resolved polygons.
+            plate_partitioner = pygplates.PlatePartitioner(
+                resolved_topologies, self.rotation_model
+            )
+        else:
+            # Some of 'curr_points' will be None so 'curr_valid_points' contains only the valid (not None)
+            # points, and 'curr_valid_points_indices' is the same length as 'curr_points' but indexes into
+            # 'curr_valid_points' so we can quickly find which point (and hence which resolved topology)
+            # in 'curr_valid_points' is associated with the a particular point in 'curr_points'.
+            curr_valid_points = []
+            curr_valid_points_indices = [None] * self.num_points
+            for point_index, curr_point in enumerate(self.curr_points):
+                if curr_point is not None:
+                    curr_valid_points_indices[point_index] = len(curr_valid_points)
+                    curr_valid_points.append(curr_point)
+            # For each valid current point find the resolved topology containing it.
+            resolved_topologies_containing_curr_valid_points = (
+                _ptt.utils.points_in_polygons.find_polygons(
+                    curr_valid_points,
+                    [
+                        resolved_topology.get_resolved_boundary()
+                        for resolved_topology in resolved_topologies
+                    ],
+                    resolved_topologies,
+                )
+            )
+
+        # Iterate over all points.
+        for point_index, curr_point in enumerate(self.curr_points):
+            if curr_point is None:
+                # Current point is not currently active - so skip it.
+                self.curr_topology_plate_ids[point_index] = None
+                self.curr_resolved_plate_boundaries[point_index] = None
+                continue
+
+            # Find the plate id of the polygon that contains 'curr_point'.
+            if ReconstructByTopologies.use_plate_partitioner:
+                curr_polygon = plate_partitioner.partition_point(curr_point)
+            else:
+                curr_polygon = resolved_topologies_containing_curr_valid_points[
+                    # Index back into 'curr_valid_points' and hence also into
+                    # 'resolved_topologies_containing_curr_valid_points'.
+                    curr_valid_points_indices[point_index]
+                ]
+            self.curr_resolved_plate_boundaries[point_index] = curr_polygon
+
+            # If the polygon is None, that means (presumably) that it fell into a crack between
+            # topologies. So it will be skipped and thrown away from future iterations.
+            if curr_polygon is None:
+                self.curr_topology_plate_ids[point_index] = None
+                continue
+
+            # Set the plate ID of resolved topology containing current point.
+            self.curr_topology_plate_ids[point_index] = (
+                curr_polygon.get_feature().get_reconstruction_plate_id()
+            )
\ No newline at end of file