diff --git a/Notebooks/10-SeafloorGrids.ipynb b/Notebooks/10-SeafloorGrids.ipynb index b74e6aac..ef6605bb 100644 --- a/Notebooks/10-SeafloorGrids.ipynb +++ b/Notebooks/10-SeafloorGrids.ipynb @@ -1,649 +1,735 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "6a6c21c7", - "metadata": {}, - "source": [ - "## 10 - Seafloor Grids\n", - "\n", - "An adaptation of [agegrid-01](https://github.com/siwill22/agegrid-0.1) written by Simon Williams, Nicky Wright and John Cannon for gridding general z-values onto seafloor basin points using GPlately.\n", - "\n", - "This notebook demonstrates how to generate seafloor age and spreading-rate grids through time. The figure below, generated by this notebook, shows seafloor age at 65 Ma. You may adjust the time parameter to generate grids at different geological times.\n", - "\n", - "![Figure: seafloor_age_65Ma.png](https://gplates.github.io/gplately/latest/sphinx/html/_images/seafloor_age_65Ma.png)\n", - "\n", - "### Citation:\n", - "Simon Williams, Nicky M. Wright, John Cannon, Nicolas Flament, R. Dietmar Müller, Reconstructing seafloor age distributions in lost ocean basins, Geoscience Frontiers, Volume 12, Issue 2, 2021, Pages 769-780, ISSN 1674-9871,\n", - "https://doi.org/10.1016/j.gsf.2020.06.004." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "db118eb4", - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "#\n", - "# Setting GPLATELY_DEBUG to a value above 100 will create the following debug files\n", - "# during seafloor gridding (that can be loaded into GPlates):\n", - "#\n", - "# - A separate debug file for each time (from 'max_time' to 'min_time) containing the seed points\n", - "# reconstructed from that time to 'min_time'.\n", - "# * The debug file with \"initial_ocean_basin\" in the filename contains reconstructions of the initial ocean basin seed points, and\n", - "# * each debug file with \"mid_ocean_ridge\" in the filename contains reconstructions of the seed points created along mid-ocean ridges\n", - "# that formed at that time specified in the filename.\n", - "# - One big debug file containing the reconstructions of ALL seed points created at ALL times.\n", - "#\n", - "# These files are located in the 'debug' sub-directory of the save directory.\n", - "#\n", - "#os.environ[\"GPLATELY_DEBUG\"] = \"200\"\n", - "\n", - "import gplately\n", - "\n", - "import pygplates\n", - "import glob\n", - "import matplotlib.pyplot as plt\n", - "import cartopy.crs as ccrs\n", - "from plate_model_manager import PlateModelManager" - ] - }, - { - "cell_type": "markdown", - "id": "9a4bddee", - "metadata": {}, - "source": [ - "### Define a rotation model, topology features and continents for the `PlateReconstruction` model\n", - "There are two ways to do this. To use local files, set `use_local_files = True`. To use `PlateModelManager`, set `use_local_files = False`." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b571bf94", - "metadata": {}, - "outputs": [], - "source": [ - "use_local_files = False" - ] - }, - { - "cell_type": "markdown", - "id": "4936baad", - "metadata": {}, - "source": [ - "#### 1) Manually pointing to files" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e5319664", - "metadata": {}, - "outputs": [], - "source": [ - "if use_local_files:\n", - " # Method 1: manually point to files\n", - " # you can download the model files at \n", - " # https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\n", - " input_directory = \"./SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4\"\n", - " \n", - " rotation_model = glob.glob(os.path.join(input_directory, '*.rot'))\n", - " static_polygons = input_directory+\"/shapes_static_polygons_Merdith_et_al.gpml\"\n", - " topology_features = [\n", - " input_directory+\"/250-0_plate_boundaries_Merdith_et_al.gpml\",\n", - " input_directory+\"/410-250_plate_boundaries_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Convergence_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Divergence_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Topologies_Merdith_et_al.gpml\",\n", - " input_directory+\"/1000-410-Transforms_Merdith_et_al.gpml\",\n", - " input_directory+\"/TopologyBuildingBlocks_Merdith_et_al.gpml\"\n", - " ]\n", - " \n", - " continents = input_directory+\"/shapes_continents.gpml\"\n", - " coastlines = input_directory+\"/shapes_coastlines_Merdith_et_al_v2.gpmlz\"\n", - " COBs = None\n", - "\n", - " from pathlib import Path\n", - " if not Path(input_directory).is_dir():\n", - " raise FileNotFoundError(f\"The input directory does not exist: {input_directory}. Ensure your files are in the correct locations or download the example input files at https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\")\n", - " " - ] - }, - { - "cell_type": "markdown", - "id": "1519da7c", - "metadata": {}, - "source": [ - "#### 2) Using `PlateModelManager`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3ebf9727", - "metadata": {}, - "outputs": [], - "source": [ - "if not use_local_files:\n", - " # Method 2: use PlateModelManager\n", - " pm_manager = PlateModelManager()\n", - " plate_model = pm_manager.get_model(\"Merdith2021\", data_dir=\"plate-model-repo\")\n", - " \n", - " rotation_model = plate_model.get_rotation_model()\n", - " topology_features = plate_model.get_topologies()\n", - " static_polygons = plate_model.get_static_polygons()\n", - " \n", - " continents = plate_model.get_layer('ContinentalPolygons') " - ] - }, - { - "cell_type": "markdown", - "id": "3560443e", - "metadata": {}, - "source": [ - "## The SeafloorGrid object\n", - "...is a collection of methods to generate seafloor grids.\n", - "\n", - "The critical input parameters are:\n", - "\n", - "### Plate model parameters\n", - "* **`plate_reconstruction`**: The gplately `PlateReconstruction` object defined in the cell below. This object is a collection of methods for calculating plate tectonic stats through geological time.\n", - "* **`plot_topologies`**: The gplately `PlotTopologies` object defined in the cell below. This object is a collection of methods for resolving geological features needed for gridding to a certain geological time.\n", - "\n", - "#### Define the `PlateReconstruction` and `PlotTopologies` objects" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "12bec53e", - "metadata": {}, - "outputs": [], - "source": [ - "model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)\n", - "\n", - "gplot = gplately.PlotTopologies(model, continents=continents)" - ] - }, - { - "cell_type": "markdown", - "id": "57d6e8b2", - "metadata": {}, - "source": [ - "### Time parameters\n", - "* **`max_time`**: The first time step to start generating age grids. This is the time step where we define certain **gridding initial conditions**, which will be explained below.\n", - "* **`min_time`**: The final step to generate age grids. This is the time when recursive reconstructions starting from `max_time` stop.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "17add2d7", - "metadata": {}, - "outputs": [], - "source": [ - "max_time = 410.\n", - "min_time = 65." - ] - }, - { - "cell_type": "markdown", - "id": "703d2315", - "metadata": {}, - "source": [ - "### Ridge resolution parameters\n", - "With each reconstruction time step, mid-ocean ridge segments emerge and spread across the ocean floor. In `gplately`, these ridges are lines, or a tessellation of infinitesimal points. The spatio-temporal resolution of these ridge points can be controlled by two parameters:\n", - "* **`ridge_time_step`**: The \"delta time\" or time increment between successive resolutions of ridges, and hence successive grids. By default this is 1Myr, so grids are produced every millionth year.\n", - "* **`ridge_sampling`**: This controls the geographical resolution (in degrees) with which ridge lines are partitioned into points. The larger this is, the larger the spacing between successive ridge points, and hence the smaller the number of ridge points at each timestep. By default this is 0.5 degrees, so points are spaced 0.5 degrees apart. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b959a7b4", - "metadata": {}, - "outputs": [], - "source": [ - "ridge_time_step = 5. # We set the time step to 5 Myr to reduce the notebook’s runtime.\n", - "ridge_sampling = 0.5" - ] - }, - { - "cell_type": "markdown", - "id": "c1c20f7d", - "metadata": {}, - "source": [ - "***\n", - "### Grid resolution parameters\n", - "The ridge resolution parameters mentioned above take care of the spatio-temporal positioning of ridge features/points. However, these ridge points need to be interpolated onto regular grids - separate parameters are needed for these grids. \n", - "\n", - "#### For the `max_time` initial condition\n", - "At `max_time`, i.e. 410Ma, the `PlateReconstruction` model has not been initialised at 411Ma and any time(s) before that to sculpt the geological history before `max_time`. In terms of the workflow, all geological history before `max_time` is unknown. Thus, the global seafloor point distriubtion, and each point's spreading rate and age must be manually defined. \n", - "\n", - "The initial point distribution is an icosahedral point mesh, made using [`stripy`](https://github.com/underworldcode/stripy/tree/294354c00dd72e085a018e69c345d9353c6fafef).\n", - "\n", - "* **`refinement_levels`**: A unitless integer that controls the number of points in the `max_time` icosahedral point mesh.\n", - " \n", - " 6 is the default. Any higher will result in a finer mesh.\n", - " \n", - " \n", - "* **`initial_ocean_mean_spreading_rate`** (in units of mm/yr or km/myr): Since the geological history before and at `max_time` is unknown, we will need to manually define the spreading rates and ages of all ocean points. This is manually set to a uniform spreading rate of 75 (mm/yr or km/myr). Each point's age is equal to its proximity to the nearest mid ocean ridge (assuming that ridge is the source of the point) divided by half this uniform spreading rate (half to account for spreading to the left and right of a ridge). \n", - "\n", - "#### For the interpolated regular grids\n", - "The regular grid on which all data is interpolated has a resolution that can be controlled by:\n", - "\n", - "* **`grid_spacing`**: The degree spacing between successive nodes in the grid. By default, this is 0.1 degrees. **Acceptable degree spacings are 0.1, 0.25, 0.5, 0.75, or 1 degree(s)** because these allow cleanly divisible grid nodes across global latitude and longitude extents. Anything greater than 1 will be rounded to 1, and anything between these acceptable spacings will be rounded to the nearest acceptable spacing. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5ec00f36", - "metadata": {}, - "outputs": [], - "source": [ - "refinement_levels = 6\n", - "initial_ocean_mean_spreading_rate = 50.\n", - "\n", - "# Gridding parameter\n", - "grid_spacing = 0.25\n", - "\n", - "extent=[-180,180,-90,90]" - ] - }, - { - "cell_type": "markdown", - "id": "8464bfb4", - "metadata": {}, - "source": [ - "***\n", - "### File saving and naming parameters\n", - "When grids are produced, they are saved to:\n", - "* **`save_directory`**: A string to a directory that must exist already.\n", - "\n", - "These files are named according to:\n", - "* **`file_collection`**: A string used to help with naming all output files. \n", - "***" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "106b61f5", - "metadata": {}, - "outputs": [], - "source": [ - "# continent masks, initial ocean seed points, and gridding input files are kept here\n", - "output_parent_directory = os.path.join(\n", - " \"NotebookFiles\",\n", - " \"Notebook10\",\n", - ")\n", - "save_directory = os.path.join(\n", - " output_parent_directory,\n", - " \"seafloor_grid_output\",\n", - ")\n", - "print(f\"The ouput files created by this notebook will be saved in folder {save_directory}.\" )\n", - "\n", - "# A string to help name files according to a plate model \"Merdith2021\"\n", - "file_collection = \"Merdith2021\"\n" - ] - }, - { - "cell_type": "markdown", - "id": "a328fa6d", - "metadata": {}, - "source": [ - "***\n", - "### Continent collision\n", - "\n", - "One way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is collides into continental crust. At each timestep continental crust is represented as a **continental mask** which is a binary grid that partitions continental regions from oceanic regions. By default, continental masks are generated from the reconstructed continental polygons (where the continent polygons features are specified with ``continent_polygon_features``). But they can also be generated using [continent contouring](https://github.com/EarthByte/continent-contouring) (by specifying ``use_continent_contouring=True``) or you can explicitly provide them as your own continent masks (by specifying them with ``continent_mask_filename``).\n", - "\n", - "### Subduction parameters\n", - "Another way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is approaching a plate boundary such that a **velocity and displacement test** is passed:\n", - "\n", - "#### Velocity test\n", - "First, given the trajectory of a point at the next time step, a point will cross a subducting plate boundary (and thus will be deleted) if the difference between its velocity on its current plate AND the velocity it will have on the other plate is greater than this velocity difference is higher than a specified **threshold delta velocity in kms/Myr**. \n", - "\n", - "#### Displacement test\n", - "If the proximity of the point's previous time position to the plate boundary it is approaching is higher than a set distance threshold (**threshold distance to boundary per Myr in kms/Myr**), then the point is far enough away from the boundary that it cannot be subducted or consumed by it, and \n", - "hence the point is still active.\n", - "\n", - "The parameter needed to encapsulate these thresholds is:\n", - "* **`subduction_collision_parameters`**: This a tuple with two elements, by default (5.0, 10.0), for the (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my)\n", - "***" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "68450bf6", - "metadata": {}, - "outputs": [], - "source": [ - "continent_polygon_features=continents\n", - "subduction_collision_parameters=(5.0, 10.0)" - ] - }, - { - "cell_type": "markdown", - "id": "bd978b09", - "metadata": {}, - "source": [ - "***\n", - "### Methodology parameters\n", - "* **`resume_from_checkpoints`**: Gridding can take a while (depending on the number of CPUs used). If this is changed to `True`, after gridding was interrupted, then rerunning the gridding cell will resume from where it left off.\n", - "* **`nprocs`**: The number of CPUs to use for parts of the code that are parallelized. The default is `-2` which uses all CPUs except one (to keep the system responsive).\n", - "***" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b1a6fbee", - "metadata": {}, - "outputs": [], - "source": [ - "# Methodology parameters\n", - "resume_from_checkpoints = False\n", - "nprocs = -2" - ] - }, - { - "cell_type": "markdown", - "id": "5038c7e6", - "metadata": {}, - "source": [ - "Use all parameters to define `SeafloorGrid`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ed5fc9e9", - "metadata": {}, - "outputs": [], - "source": [ - "# The SeafloorGrid object with all aforementioned parameters\n", - "seafloorgrid = gplately.SeafloorGrid(\n", - " \n", - " plate_reconstruction = model, \n", - " \n", - " # Time parameters\n", - " max_time = max_time,\n", - " min_time = min_time,\n", - " \n", - " # Ridge tessellation parameters\n", - " ridge_time_step = ridge_time_step,\n", - " ridge_sampling = ridge_sampling,\n", - " \n", - " # Gridding parameters\n", - " grid_spacing = grid_spacing,\n", - "\n", - " extent = extent,\n", - " \n", - " # Naming parameters\n", - " save_directory = save_directory,\n", - " file_collection = file_collection,\n", - " \n", - " # Initial condition parameters\n", - " refinement_levels = refinement_levels,\n", - " initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate,\n", - " \n", - " # Subduction parameters\n", - " subduction_collision_parameters = subduction_collision_parameters,\n", - " \n", - " # Methodology parameters\n", - " resume_from_checkpoints = resume_from_checkpoints,\n", - "\n", - " # Continental polygons.\n", - " continent_polygon_features = continent_polygon_features,\n", - "\n", - " nprocs = nprocs\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "ab3cbf87", - "metadata": {}, - "source": [ - "## How to use `SeafloorGrid`\n", - "\n", - "### 1) Run `seafloorgrid.reconstruct_by_topologies()`\n", - "Running `seafloorgrid.reconstruct_by_topologies()` prepares all ocean basin points and their seafloor ages and spreading rates for gridding, and reconstructs all active points per timestep to form the grids per timestep.\n", - "\n", - "#### Continental masks\n", - "A netCDF4 continent mask is produced for each _time_ from `max_time` to `min_time`.\n", - "\n", - "#### Initial ocean basin\n", - "At `max_time`, an initial ocean seed point icosahedral mesh fills the ocean basin. Each initial seed point is allocated an age based on its distance from the nearest mid-ocean ridge. And these are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at `max_time`).\n", - "\n", - "#### Mid-ocean ridges\n", - "At each successive timestep (`ridge_time_step`), new seed points emerge from spreading ridges, and they are allocated their own spreading rates. The seed points emerging at a particular _time_ are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at that _time_). This results in one `.npz` file for each _time_ from `max_time - ridge_time_step` to `min_time`.\n", - "\n", - "#### Seafloor age and spreading rate\n", - "After all seed points (initial and mid-ocean ridge) have been reconstructed to `min_time`, then one `.npz` file is generated for each _time_ from `max_time` to `min_time`. These files contain the gridding input that will later be used to generate the seafloor age and spreading rate grids.\n", - "\n", - "#### Resuming an interrupted gridding run\n", - "*If the `resume_from_checkpoints` parameter is passed as `True`, after `seafloorgrid.reconstruct_by_topologies()` was interrupted, then you can rerun the cell. The workflow will pick up from where it left off (provided the continent mask files and reconstructed seed point data files have not been erased).*\n", - "\n", - "To overwrite all files in `save_directory`, pass `resume_from_checkpoints` as `False` (this is the default behaviour).\n", - "\n", - "#### Reconstruct by topologies\n", - "The `ReconstructByTopologies` object (written by Simon Williams, Nicky Wright, and John Cannon) handles the reconstruction of seed points. `ReconstructByTopologies` (RBT) identifies active points on the ocean basin per timestep. It works as follows:\n", - "\n", - "If an ocean point on one plate ID transitions into another rigid plate ID at the next timestep, RBT calculates the point's velocity difference between both plates. The point **may** have subducted/collided with a continent at this boundary if this velocity difference is higher than a set velocity threshold. To ascertain whether the point should indeed be deactivated, a second test is conducted: RBT checks the previous time position of the point and calculates this point’s proximity to the boundary of the plate ID polygon it is approaching. If this distance is higher than a set distance 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. Else, it is deactivated/deleted.\n", - "\n", - "Once all active points and their z-values are identified, they are written to the gridding input file (.npz) for that timestep.\n", - "\n", - "Below is an example of the initial condition: the ocean basin is populated with an `intial_mean_ocean_spreading_rate` of 50 mm/yr at a `max_time` of 410Ma. Reconstruction over 10 Myr to 400 Ma sees new seed points emerging from ridges with their own plate-model-ascribed spreading rates.\n", - "\n", - "![init_condition](https://raw.githubusercontent.com/GPlates/gplately/master/Notebooks/NotebookFiles/Notebook10/Merdith2021_initial_sr_conditions.png)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "143ea077", - "metadata": {}, - "outputs": [], - "source": [ - "import datetime\n", - "import time\n", - "\n", - "start = time.time()\n", - "seafloorgrid.reconstruct_by_topologies()\n", - "end = time.time()\n", - "duration = datetime.timedelta(seconds=end - start)\n", - "print(\"Duration: {}\".format(duration))" - ] - }, - { - "cell_type": "markdown", - "id": "a0a0ca60", - "metadata": {}, - "source": [ - "### 2) Run `seafloorgrid.lat_lon_z_to_netCDF` to write grids to netCDF\n", - "Calling `seafloorgrid.lat_lon_z_to_netCDF` grids one set of z-data per latitude-longitude pair from each timestep's gridding input file (produced in `seafloorgrid.reconstruct_by_topologies()`). Grids are in netCDF format.\n", - "\n", - "The desired z-data to grid is identified using a `zval_name`. \n", - "For example, seafloor age grids can be produced using `SEAFLOOR_AGE`, and spreading rate grids are `SPREADING_RATE`. \n", - "\n", - "Use `unmasked = True` to output both the masked and unmasked versions of the grids." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "891183f3", - "metadata": {}, - "outputs": [], - "source": [ - "seafloorgrid.lat_lon_z_to_netCDF(\"SEAFLOOR_AGE\", unmasked=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d51a41a5", - "metadata": {}, - "outputs": [], - "source": [ - "seafloorgrid.lat_lon_z_to_netCDF(\"SPREADING_RATE\", unmasked=True)" - ] - }, - { - "cell_type": "markdown", - "id": "8213b46a", - "metadata": {}, - "source": [ - "### Plotting a sample age grid and spreading rate grid\n", - "Read one netCDF grid using GPlately's `Raster` object from `grids`, and plot it using the `PlotTopologies` object.\n", - "\n", - "Notice the evolution of seafloor spreading rate from the initial value set with `initial_ocean_mean_spreading_rate`. Eventually, this initial uniform spreading rate will be phased out with sufficient recursive reconstruction." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "06a92ce9-ecea-49b9-9c99-c33438e75379", - "metadata": {}, - "outputs": [], - "source": [ - "time = min_time\n", - "# set the plot_unmasked_grid variable to True if you would like to see the unmasked age grid.\n", - "plot_unmasked_grid = False\n", - "\n", - "# age grids\n", - "if not plot_unmasked_grid:\n", - " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_{time:0.2f}Ma.nc\"\n", - "else:\n", - " agegrid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_unmasked_{time:0.2f}Ma.nc\"\n", - "\n", - "age_grid = gplately.grids.Raster(agegrid_filename)\n", - "\n", - "# Prepare plots\n", - "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", - "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", - "\n", - "gplot.time = time\n", - "plt.title(f\"Seafloor Age at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", - "\n", - "agegrid_cmap = \"YlGnBu\"\n", - "agegrid_cpt_file = \"NotebookFiles/agegrid.cpt\"\n", - "if not os.path.isfile(agegrid_cpt_file):\n", - " print(f\"The file {agegrid_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/agegrid.cpt\")\n", - " print(f\"The '{agegrid_cmap}' will be used to plot the age grid.\")\n", - " agegrid_cpt_file = None\n", - "if agegrid_cpt_file:\n", - " try:\n", - " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", - " agegrid_cmap=get_cmap_from_gmt_cpt(agegrid_cpt_file)\n", - " except ImportError as e:\n", - " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", - " print(e)\n", - "\n", - "im = gplot.plot_grid(\n", - " ax, \n", - " age_grid.data, \n", - " cmap = agegrid_cmap,\n", - " vmin = 0, \n", - " vmax =410,\n", - ")\n", - "gplot.plot_ridges(ax,linewidth=1)\n", - "if not plot_unmasked_grid:\n", - " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", - "plt.colorbar(im, label=\"Seafloor Age (Myr)\", shrink=0.5, pad=0.05)\n", - "\n", - "# Save figure\n", - "plt.savefig(\n", - " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\",\n", - " bbox_inches = 'tight',\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "56e03e39", - "metadata": {}, - "outputs": [], - "source": [ - "time = min_time\n", - "# set the plot_unmasked_grid variable to True if you would like to see the unmasked spreading rate grid.\n", - "plot_unmasked_grid = False\n", - "\n", - "# spreading rate grids\n", - "if not plot_unmasked_grid:\n", - " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_{time:0.2f}Ma.nc\"\n", - "else:\n", - " srgrid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_unmasked_{time:0.2f}Ma.nc\"\n", - "\n", - "sr_grid = gplately.grids.Raster(srgrid_filename)\n", - "\n", - "# Prepare plots\n", - "fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", - "ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", - "\n", - "gplot.time = time\n", - "plt.title(f\"Seafloor Spreading Rate at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", - "\n", - "spreading_rate_cmap = \"magma\"\n", - "spreading_rate_cpt_file = \"NotebookFiles/spreading_full_rate.cpt\"\n", - "if not os.path.isfile(spreading_rate_cpt_file):\n", - " print(f\"The file {spreading_rate_cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/NotebookFiles/spreading_full_rate.cpt\")\n", - " print(f\"The '{spreading_rate_cmap}' will be used to plot the spreading rate grid.\")\n", - " spreading_rate_cpt_file = None\n", - "if spreading_rate_cpt_file:\n", - " try:\n", - " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", - " spreading_rate_cmap=get_cmap_from_gmt_cpt(spreading_rate_cpt_file)\n", - " except ImportError as e:\n", - " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", - " print(e)\n", - "\n", - "im = gplot.plot_grid(\n", - " ax, \n", - " sr_grid.data, \n", - " cmap = spreading_rate_cmap,\n", - " vmin = 0, \n", - " vmax =200,\n", - ")\n", - "gplot.plot_ridges(ax, linewidth=1)\n", - "if not plot_unmasked_grid:\n", - " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", - "plt.colorbar(im, label=\"Spreading rate (mm/yr)\", shrink=0.5, pad=0.05, ticks=[0,50,100,150,200])\n", - "\n", - "# Save figure\n", - "plt.savefig(\n", - " f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\",\n", - " bbox_inches = 'tight',\n", - ")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.11.11" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6a6c21c7", + "metadata": {}, + "source": [ + "## 10 - Seafloor Grids\n", + "\n", + "An adaptation of [agegrid-01](https://github.com/siwill22/agegrid-0.1) written by Simon Williams, Nicky Wright and John Cannon for gridding general z-values onto seafloor basin points using GPlately.\n", + "\n", + "This notebook demonstrates how to generate seafloor age and spreading-rate grids through time. The figure below, generated by this notebook, shows seafloor age at 65 Ma. You may adjust the time parameter to generate grids at different geological times.\n", + "\n", + "![Figure: seafloor_age_65Ma.png](https://gplates.github.io/gplately/latest/sphinx/html/_images/seafloor_age_65Ma.png)\n", + "\n", + "### Citation:\n", + "Simon Williams, Nicky M. Wright, John Cannon, Nicolas Flament, R. Dietmar Müller, Reconstructing seafloor age distributions in lost ocean basins, Geoscience Frontiers, Volume 12, Issue 2, 2021, Pages 769-780, ISSN 1674-9871,\n", + "https://doi.org/10.1016/j.gsf.2020.06.004." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db118eb4", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "#\n", + "# Setting GPLATELY_DEBUG to a value above 100 will create the following debug files\n", + "# during seafloor gridding (that can be loaded into GPlates):\n", + "#\n", + "# - A separate debug file for each time (from 'max_time' to 'min_time) containing the seed points\n", + "# reconstructed from that time to 'min_time'.\n", + "# * The debug file with \"initial_ocean_basin\" in the filename contains reconstructions of the initial ocean basin seed points, and\n", + "# * each debug file with \"mid_ocean_ridge\" in the filename contains reconstructions of the seed points created along mid-ocean ridges\n", + "# that formed at that time specified in the filename.\n", + "# - One big debug file containing the reconstructions of ALL seed points created at ALL times.\n", + "#\n", + "# These files are located in the 'debug' sub-directory of the save directory.\n", + "#\n", + "#os.environ[\"GPLATELY_DEBUG\"] = \"200\"\n", + "\n", + "import gplately\n", + "\n", + "import glob\n", + "import matplotlib.pyplot as plt\n", + "import cartopy.crs as ccrs\n", + "from plate_model_manager import PlateModelManager" + ] + }, + { + "cell_type": "markdown", + "id": "9a4bddee", + "metadata": {}, + "source": [ + "### Define a rotation model, topology features and continents for the `PlateReconstruction` model\n", + "There are two ways to do this. To use local files, set `use_local_files = True`. To use `PlateModelManager`, set `use_local_files = False`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b571bf94", + "metadata": {}, + "outputs": [], + "source": [ + "use_local_files = False" + ] + }, + { + "cell_type": "markdown", + "id": "4936baad", + "metadata": {}, + "source": [ + "#### 1) Manually pointing to files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5319664", + "metadata": {}, + "outputs": [], + "source": [ + "if use_local_files:\n", + " # Method 1: manually point to files\n", + " # you can download the model files at \n", + " # https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\n", + " input_directory = \"./SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4\"\n", + " \n", + " rotation_model = glob.glob(os.path.join(input_directory, '*.rot'))\n", + " static_polygons = input_directory+\"/shapes_static_polygons_Merdith_et_al.gpml\"\n", + " topology_features = [\n", + " input_directory+\"/250-0_plate_boundaries_Merdith_et_al.gpml\",\n", + " input_directory+\"/410-250_plate_boundaries_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Convergence_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Divergence_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Topologies_Merdith_et_al.gpml\",\n", + " input_directory+\"/1000-410-Transforms_Merdith_et_al.gpml\",\n", + " input_directory+\"/TopologyBuildingBlocks_Merdith_et_al.gpml\"\n", + " ]\n", + " \n", + " continents = input_directory+\"/shapes_continents.gpml\"\n", + " coastlines = input_directory+\"/shapes_coastlines_Merdith_et_al_v2.gpmlz\"\n", + " COBs = None\n", + "\n", + " from pathlib import Path\n", + " if not Path(input_directory).is_dir():\n", + " raise FileNotFoundError(f\"The input directory does not exist: {input_directory}. Ensure your files are in the correct locations or download the example input files at https://www.earthbyte.org/webdav/ftp/Data_Collections/Merdith_etal_2021_ESR/SM2-Merdith_et_al_1_Ga_reconstruction_v1.2.4.zip\")\n", + " " + ] + }, + { + "cell_type": "markdown", + "id": "1519da7c", + "metadata": {}, + "source": [ + "#### 2) Using `PlateModelManager`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3ebf9727", + "metadata": {}, + "outputs": [], + "source": [ + "if not use_local_files:\n", + " # Method 2: use PlateModelManager\n", + " pm_manager = PlateModelManager()\n", + " plate_model = pm_manager.get_model(\"Merdith2021\", data_dir=\"plate-model-repo\")\n", + " \n", + " rotation_model = plate_model.get_rotation_model()\n", + " topology_features = plate_model.get_topologies()\n", + " static_polygons = plate_model.get_static_polygons()\n", + " \n", + " continents = plate_model.get_layer('ContinentalPolygons') " + ] + }, + { + "cell_type": "markdown", + "id": "3560443e", + "metadata": {}, + "source": [ + "## The SeafloorGrid object\n", + "...is a collection of methods to generate seafloor grids.\n", + "\n", + "The critical input parameters are:\n", + "\n", + "### Plate model parameters\n", + "* **`plate_reconstruction`**: The gplately `PlateReconstruction` object defined in the cell below. This object is a collection of methods for calculating plate tectonic stats through geological time.\n", + "* **`plot_topologies`**: The gplately `PlotTopologies` object defined in the cell below. This object is a collection of methods for resolving geological features needed for gridding to a certain geological time.\n", + "\n", + "#### Define the `PlateReconstruction` and `PlotTopologies` objects" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12bec53e", + "metadata": {}, + "outputs": [], + "source": [ + "model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)\n", + "\n", + "gplot = gplately.PlotTopologies(model, continents=continents)" + ] + }, + { + "cell_type": "markdown", + "id": "57d6e8b2", + "metadata": {}, + "source": [ + "### Time parameters\n", + "* **`max_time`**: The first time step to start generating age grids. This is the time step where we define certain **gridding initial conditions**, which will be explained below.\n", + "* **`min_time`**: The final step to generate age grids. This is the time when recursive reconstructions starting from `max_time` stop.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17add2d7", + "metadata": {}, + "outputs": [], + "source": [ + "max_time = 410.\n", + "min_time = 65." + ] + }, + { + "cell_type": "markdown", + "id": "703d2315", + "metadata": {}, + "source": [ + "### Ridge resolution parameters\n", + "With each reconstruction time step, mid-ocean ridge segments emerge and spread across the ocean floor. In `gplately`, these ridges are lines, or a tessellation of infinitesimal points. The spatio-temporal resolution of these ridge points can be controlled by two parameters:\n", + "* **`ridge_time_step`**: The \"delta time\" or time increment between successive resolutions of ridges, and hence successive grids. By default this is 1Myr, so grids are produced every millionth year.\n", + "* **`ridge_sampling`**: This controls the geographical resolution (in degrees) with which ridge lines are partitioned into points. The larger this is, the larger the spacing between successive ridge points, and hence the smaller the number of ridge points at each timestep. By default this is 0.5 degrees, so points are spaced 0.5 degrees apart. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b959a7b4", + "metadata": {}, + "outputs": [], + "source": [ + "ridge_time_step = 5. # We set the time step to 5 Myr to reduce the notebook’s runtime.\n", + "ridge_sampling = 0.5" + ] + }, + { + "cell_type": "markdown", + "id": "c1c20f7d", + "metadata": {}, + "source": [ + "***\n", + "### Grid resolution parameters\n", + "The ridge resolution parameters mentioned above take care of the spatio-temporal positioning of ridge features/points. However, these ridge points need to be interpolated onto regular grids - separate parameters are needed for these grids. \n", + "\n", + "#### For the `max_time` initial condition\n", + "At `max_time`, i.e. 410Ma, the `PlateReconstruction` model has not been initialised at 411Ma and any time(s) before that to sculpt the geological history before `max_time`. In terms of the workflow, all geological history before `max_time` is unknown. Thus, the global seafloor point distriubtion, and each point's spreading rate and age must be manually defined. \n", + "\n", + "The initial point distribution is an icosahedral point mesh, made using [`stripy`](https://github.com/underworldcode/stripy/tree/294354c00dd72e085a018e69c345d9353c6fafef).\n", + "\n", + "* **`refinement_levels`**: A unitless integer that controls the number of points in the `max_time` icosahedral point mesh.\n", + " \n", + " 6 is the default. Any higher will result in a finer mesh.\n", + " \n", + " \n", + "* **`initial_ocean_mean_spreading_rate`** (in units of mm/yr or km/myr): Since the geological history before and at `max_time` is unknown, we will need to manually define the spreading rates and ages of all ocean points. This is manually set to a uniform spreading rate of 75 (mm/yr or km/myr). Each point's age is equal to its proximity to the nearest mid ocean ridge (assuming that ridge is the source of the point) divided by half this uniform spreading rate (half to account for spreading to the left and right of a ridge). \n", + "\n", + "#### For the interpolated regular grids\n", + "The regular grid on which all data is interpolated has a resolution that can be controlled by:\n", + "\n", + "* **`grid_spacing`**: The degree spacing between successive nodes in the grid. By default, this is 0.1 degrees. **Acceptable degree spacings are 0.1, 0.25, 0.5, 0.75, or 1 degree(s)** because these allow cleanly divisible grid nodes across global latitude and longitude extents. Anything greater than 1 will be rounded to 1, and anything between these acceptable spacings will be rounded to the nearest acceptable spacing. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ec00f36", + "metadata": {}, + "outputs": [], + "source": [ + "refinement_levels = 6\n", + "initial_ocean_mean_spreading_rate = 75.\n", + "\n", + "# Gridding parameter\n", + "grid_spacing = 0.25\n", + "\n", + "extent=[-180,180,-90,90]" + ] + }, + { + "cell_type": "markdown", + "id": "8464bfb4", + "metadata": {}, + "source": [ + "***\n", + "### File saving and naming parameters\n", + "When grids are produced, they are saved to:\n", + "* **`save_directory`**: A string to a directory that must exist already.\n", + "\n", + "These files are named according to:\n", + "* **`file_collection`**: A string used to help with naming all output files. \n", + "***" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "106b61f5", + "metadata": {}, + "outputs": [], + "source": [ + "# continent masks, initial ocean seed points, and gridding input files are kept here\n", + "output_parent_directory = os.path.join(\n", + " \"NotebookFiles\",\n", + " \"Notebook10\",\n", + ")\n", + "save_directory = os.path.join(\n", + " output_parent_directory,\n", + " \"seafloor_grid_output\",\n", + ")\n", + "print(f\"The ouput files created by this notebook will be saved in folder {save_directory}.\" )\n", + "\n", + "# A string to help name files according to a plate model \"Merdith2021\"\n", + "file_collection = \"Merdith2021\"\n" + ] + }, + { + "cell_type": "markdown", + "id": "a328fa6d", + "metadata": {}, + "source": [ + "***\n", + "### Continent collision\n", + "\n", + "One way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is collides into continental crust. At each timestep continental crust is represented as a **continental mask** which is a binary grid that partitions continental regions from oceanic regions. By default, continental masks are generated from the reconstructed continental polygons (where the continent polygons features are specified with ``continent_polygon_features``). But they can also be generated using [continent contouring](https://github.com/EarthByte/continent-contouring) (by specifying ``use_continent_contouring=True``) or you can explicitly provide them as your own continent masks (by specifying them with ``continent_mask_filename``).\n", + "\n", + "### Subduction collision\n", + "Another way an oceanic point can be deleted from the ocean basin at a certain timestep is if it is approaching a convergent plate boundary. This case is now handled automatically without specifying any parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68450bf6", + "metadata": {}, + "outputs": [], + "source": [ + "continent_polygon_features=continents" + ] + }, + { + "cell_type": "markdown", + "id": "bd978b09", + "metadata": {}, + "source": [ + "***\n", + "### Methodology parameters\n", + "* **`resume_from_checkpoints`**: Gridding can take a while (depending on the number of CPUs used). If this is changed to `True`, after gridding was interrupted, then rerunning the gridding cell will resume from where it left off.\n", + "* **`nprocs`**: The number of CPUs to use for parts of the code that are parallelized. The default is `-2` which uses all CPUs except one (to keep the system responsive).\n", + "***" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b1a6fbee", + "metadata": {}, + "outputs": [], + "source": [ + "# Methodology parameters\n", + "resume_from_checkpoints = False\n", + "nprocs = -2" + ] + }, + { + "cell_type": "markdown", + "id": "5038c7e6", + "metadata": {}, + "source": [ + "Use all parameters to define `SeafloorGrid`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed5fc9e9", + "metadata": {}, + "outputs": [], + "source": [ + "# The SeafloorGrid object with all aforementioned parameters\n", + "seafloorgrid = gplately.SeafloorGrid(\n", + " \n", + " model,\n", + " \n", + " # Time parameters\n", + " max_time = max_time,\n", + " min_time = min_time,\n", + " \n", + " # Ridge tessellation parameters\n", + " ridge_time_step = ridge_time_step,\n", + " ridge_sampling = ridge_sampling,\n", + " \n", + " # Gridding parameters\n", + " grid_spacing = grid_spacing,\n", + "\n", + " extent = extent,\n", + " \n", + " # Naming parameters\n", + " save_directory = save_directory,\n", + " file_collection = file_collection,\n", + " \n", + " # Initial condition parameters\n", + " refinement_levels = refinement_levels,\n", + " initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate,\n", + " \n", + " # Methodology parameters\n", + " resume_from_checkpoints = resume_from_checkpoints,\n", + "\n", + " # Continental polygons.\n", + " continent_polygon_features = continent_polygon_features,\n", + "\n", + " nprocs = nprocs\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "ab3cbf87", + "metadata": {}, + "source": [ + "## How to use `SeafloorGrid`\n", + "\n", + "### 1) Run `seafloorgrid.reconstruct_by_topologies()`\n", + "Running `seafloorgrid.reconstruct_by_topologies()` prepares all ocean basin points and their seafloor ages and spreading rates for gridding, and reconstructs all active points per timestep to form the grids per timestep.\n", + "\n", + "#### Continental masks\n", + "A netCDF4 continent mask is produced for each _time_ from `max_time` to `min_time`.\n", + "\n", + "#### Initial ocean basin\n", + "At `max_time`, an initial ocean seed point icosahedral mesh fills the ocean basin. Each initial seed point is allocated an age based on its distance from the nearest mid-ocean ridge. And these are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at `max_time`).\n", + "\n", + "#### Mid-ocean ridges\n", + "At each successive timestep (`ridge_time_step`), new seed points emerge from spreading ridges, and they are allocated their own spreading rates. The seed points emerging at a particular _time_ are reconstructed to `min_time` to produce an `.npz` file (containing reconstructed seed point data for seed points created at that _time_). This results in one `.npz` file for each _time_ from `max_time - ridge_time_step` to `min_time`.\n", + "\n", + "#### Seafloor age and spreading rate\n", + "After all seed points (initial and mid-ocean ridge) have been reconstructed to `min_time`, then one `.npz` file is generated for each _time_ from `max_time` to `min_time`. These files contain the gridding input that will later be used to generate the seafloor age and spreading rate grids.\n", + "\n", + "#### Resuming an interrupted gridding run\n", + "*If the `resume_from_checkpoints` parameter is passed as `True`, after `seafloorgrid.reconstruct_by_topologies()` was interrupted, then you can rerun the cell. The workflow will pick up from where it left off (provided the continent mask files and reconstructed seed point data files have not been erased).*\n", + "\n", + "To overwrite all files in `save_directory`, pass `resume_from_checkpoints` as `False` (this is the default behaviour).\n", + "\n", + "#### Reconstruct by topologies\n", + "The `ReconstructByTopologies` object (written by Simon Williams, Nicky Wright, and John Cannon) handles the reconstruction of seed points. `ReconstructByTopologies` (RBT) identifies active points on the ocean basin per timestep. It works as follows:\n", + "\n", + "If an ocean point on one plate ID transitions into another rigid plate ID at the next timestep, RBT calculates the point's velocity difference between both plates. The point **may** have subducted/collided with a continent at this boundary if this velocity difference is higher than a set velocity threshold. To ascertain whether the point should indeed be deactivated, a second test is conducted: RBT checks the previous time position of the point and calculates this point’s proximity to the boundary of the plate ID polygon it is approaching. If this distance is higher than a set distance 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. Else, it is deactivated/deleted.\n", + "\n", + "Once all active points and their z-values are identified, they are written to the gridding input file (.npz) for that timestep.\n", + "\n", + "Below is an example of the initial condition: the ocean basin is populated with an `intial_mean_ocean_spreading_rate` of 50 mm/yr at a `max_time` of 410Ma. Reconstruction over 10 Myr to 400 Ma sees new seed points emerging from ridges with their own plate-model-ascribed spreading rates.\n", + "\n", + "![init_condition](https://raw.githubusercontent.com/GPlates/gplately/master/Notebooks/NotebookFiles/Notebook10/Merdith2021_initial_sr_conditions.png)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "143ea077", + "metadata": {}, + "outputs": [], + "source": [ + "import datetime\n", + "import time\n", + "\n", + "start = time.time()\n", + "seafloorgrid.reconstruct_by_topologies()\n", + "end = time.time()\n", + "duration = datetime.timedelta(seconds=end - start)\n", + "print(\"Duration: {}\".format(duration))" + ] + }, + { + "cell_type": "markdown", + "id": "a0a0ca60", + "metadata": {}, + "source": [ + "### 2) Run `seafloorgrid.lat_lon_z_to_netCDF` to write grids to netCDF\n", + "Calling `seafloorgrid.lat_lon_z_to_netCDF` grids one set of z-data per latitude-longitude pair from each timestep's gridding input file (produced in `seafloorgrid.reconstruct_by_topologies()`). Grids are in netCDF format.\n", + "\n", + "The desired z-data to grid is identified using a `zval_name`. \n", + "For example, seafloor age grids can be produced using `SEAFLOOR_AGE`, and spreading rate grids are `SPREADING_RATE`. \n", + "\n", + "Use `unmasked = True` to output both the masked and unmasked versions of the grids." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "891183f3", + "metadata": {}, + "outputs": [], + "source": [ + "seafloorgrid.lat_lon_z_to_netCDF(\"SEAFLOOR_AGE\", unmasked=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d51a41a5", + "metadata": {}, + "outputs": [], + "source": [ + "seafloorgrid.lat_lon_z_to_netCDF(\"SPREADING_RATE\", unmasked=True)" + ] + }, + { + "cell_type": "markdown", + "id": "8213b46a", + "metadata": {}, + "source": [ + "### Plotting a sample age grid and spreading rate grid\n", + "Read one netCDF grid using GPlately's `Raster` object from `grids`, and plot it using the `PlotTopologies` object.\n", + "\n", + "Notice the evolution of seafloor spreading rate from the initial value set with `initial_ocean_mean_spreading_rate`. Eventually, this initial uniform spreading rate will be phased out with sufficient recursive reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06a92ce9-ecea-49b9-9c99-c33438e75379", + "metadata": {}, + "outputs": [], + "source": [ + "def plot(time, zval_name, save_fig=False):\n", + " \n", + " if zval_name != \"SEAFLOOR_AGE\" and zval_name != \"SPREADING_RATE\":\n", + " raise ValueError(f\"Unknown zval_name '{zval_name}'\")\n", + " \n", + " # set the plot_unmasked_grid variable to True if you would like to see the unmasked age/spreading-rate grid.\n", + " plot_unmasked_grid = False\n", + " \n", + " # age/spreading-rate grids\n", + " if not plot_unmasked_grid:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_{time:0.2f}Ma.nc\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_{time:0.2f}Ma.nc\"\n", + " else:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SEAFLOOR_AGE/{seafloorgrid.file_collection}_SEAFLOOR_AGE_grid_unmasked_{time:0.2f}Ma.nc\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " grid_filename = f\"{seafloorgrid.save_directory}/SPREADING_RATE/{seafloorgrid.file_collection}_SPREADING_RATE_grid_unmasked_{time:0.2f}Ma.nc\"\n", + " \n", + " grid = gplately.grids.Raster(grid_filename)\n", + " \n", + " # Prepare plots\n", + " fig = plt.figure(figsize=(8,4), dpi=240, linewidth=1)\n", + " ax = fig.add_subplot(111, projection=ccrs.Mollweide(central_longitude=20))\n", + " \n", + " gplot.time = time\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " plt.title(f\"Seafloor Age at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", + " cmap = \"YlGnBu\"\n", + " cpt_file = \"NotebookFiles/agegrid.cpt\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " plt.title(f\"Seafloor Spreading Rate at {int(time)} Ma (Model:{seafloorgrid.file_collection})\")\n", + " cmap = \"magma\"\n", + " cpt_file = \"NotebookFiles/spreading_full_rate.cpt\"\n", + " \n", + " if not os.path.isfile(cpt_file):\n", + " print(f\"The file {cpt_file} does not exist locally. It can be found at https://raw.githubusercontent.com/GPlates/gplately/refs/heads/master/Notebooks/{cpt_file}\")\n", + " print(f\"The '{cmap}' will be used to plot the grid.\")\n", + " cpt_file = None\n", + " if cpt_file:\n", + " try:\n", + " from gplately.mapping.gmt_cpt import get_cmap_from_gmt_cpt\n", + " cmap=get_cmap_from_gmt_cpt(cpt_file)\n", + " except ImportError as e:\n", + " print(\"To use GMT cpt file, ensure you are using GPlately version 2.1 or higher.\")\n", + " print(e)\n", + "\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " vmin, vmax = 0, 410\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " vmin, vmax = 0, 200\n", + " \n", + " im = gplot.plot_grid(\n", + " ax, \n", + " grid.data, \n", + " cmap = cmap,\n", + " vmin = vmin, \n", + " vmax = vmax,\n", + " )\n", + " gplot.plot_ridges(ax,linewidth=1)\n", + " if not plot_unmasked_grid:\n", + " gplot.plot_continents(ax, edgecolor=\"grey\", facecolor=\"lightgrey\", linewidth=0.5)\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " plt.colorbar(im, label=\"Seafloor Age (Myr)\", shrink=0.5, pad=0.05, ticks=[0,100,200,300,400])\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " plt.colorbar(im, label=\"Spreading rate (mm/yr)\", shrink=0.5, pad=0.05, ticks=[0,50,100,150,200])\n", + " \n", + " if save_fig:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\"\n", + " \n", + " # Save figure.\n", + " plt.savefig(fig_filename, bbox_inches = 'tight')\n", + " else:\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87cd68e6-fb41-42e8-8d3a-4ee317ed002d", + "metadata": {}, + "outputs": [], + "source": [ + "plot(min_time, zval_name=\"SEAFLOOR_AGE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08979858-bd12-4cf5-accd-b372405ecc21", + "metadata": {}, + "outputs": [], + "source": [ + "plot(min_time, zval_name=\"SPREADING_RATE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38d323c9-772e-4151-9f0d-66065b6f4a73", + "metadata": {}, + "outputs": [], + "source": [ + "def save_figures(\n", + " zval_name,\n", + " min_time=min_time,\n", + " max_time=max_time,\n", + " time_step=ridge_time_step):\n", + " \n", + " from joblib import Parallel, delayed\n", + " import numpy as np\n", + "\n", + " time_range = np.arange(max_time, min_time-0.5*time_step, -time_step)\n", + " \n", + " use_parallel = True\n", + " if use_parallel:\n", + " # Produce plots in a parallel routine - a progress bar and time taken will be shown.\n", + " parallel_plots = Parallel(n_jobs=nprocs, verbose=1)(\n", + " delayed(plot) (time, zval_name, save_fig=True) for time in time_range)\n", + " \n", + " else:\n", + " for time in time_range:\n", + " plot(time, save_fig=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eb5be45a-2725-47cd-be7c-f512a25574ca", + "metadata": {}, + "outputs": [], + "source": [ + "def save_video(\n", + " zval_name,\n", + " min_time=min_time,\n", + " max_time=max_time,\n", + " time_step=ridge_time_step):\n", + " \n", + " try:\n", + " import moviepy.editor as mpy # moviepy 1.x\n", + " except ImportError:\n", + " import moviepy as mpy # moviepy 2.x\n", + " import numpy as np\n", + "\n", + " if zval_name != \"SEAFLOOR_AGE\" and zval_name != \"SPREADING_RATE\":\n", + " raise ValueError(f\"Unknown zval_name '{zval_name}'\")\n", + " \n", + " time_range = np.arange(max_time, min_time-0.5*time_step, -time_step)\n", + " \n", + " frame_list = []\n", + " for time in time_range:\n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age_{int(time)}Ma.png\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " fig_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate_{int(time)}Ma.png\"\n", + " \n", + " frame_list.append(fig_filename)\n", + " \n", + " clip = mpy.ImageSequenceClip(frame_list, fps=20)\n", + " \n", + " if zval_name == \"SEAFLOOR_AGE\":\n", + " video_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_age.mp4\"\n", + " elif zval_name == \"SPREADING_RATE\":\n", + " video_filename = f\"{output_parent_directory}/{seafloorgrid.file_collection}_seafloor_spreading_rate.mp4\"\n", + " \n", + " clip.write_videofile(\n", + " video_filename,\n", + " fps=20,\n", + " codec=\"libx264\",\n", + " bitrate=\"8000k\",\n", + " audio=False,\n", + " logger=None,\n", + " ffmpeg_params=[\n", + " \"-vf\",\n", + " \"pad=ceil(iw/2)*2:ceil(ih/2)*2\",\n", + " \"-pix_fmt\",\n", + " \"yuv420p\",\n", + " ],\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e585b0cf-cd2f-4ead-b19b-4f926ea3139c", + "metadata": {}, + "outputs": [], + "source": [ + "generate_seafloor_age_video = False\n", + "\n", + "if generate_seafloor_age_video:\n", + " save_figures(zval_name=\"SEAFLOOR_AGE\")\n", + " save_video(zval_name=\"SEAFLOOR_AGE\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3202b8cc-d134-46ea-b1db-4fc747adaec8", + "metadata": {}, + "outputs": [], + "source": [ + "generate_spreading_rate_video = False\n", + "\n", + "if generate_spreading_rate_video:\n", + " save_figures(zval_name=\"SPREADING_RATE\")\n", + " save_video(zval_name=\"SPREADING_RATE\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb b/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb index ff30d24a..071da249 100644 --- a/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb +++ b/Notebooks/14-RuleBasedGPMLProcessingPipeline.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "3d561d78", + "id": "b2b60c22", "metadata": {}, "source": [ "\n", @@ -19,7 +19,7 @@ }, { "cell_type": "markdown", - "id": "4ca19dbe", + "id": "c24098e8", "metadata": {}, "source": [ "\n", @@ -30,7 +30,7 @@ }, { "cell_type": "markdown", - "id": "743d5977", + "id": "2938cc71", "metadata": { "lines_to_next_cell": 0 }, @@ -41,7 +41,7 @@ { "cell_type": "code", "execution_count": null, - "id": "680ba013", + "id": "40352ed2", "metadata": {}, "outputs": [], "source": [ @@ -118,7 +118,7 @@ }, { "cell_type": "markdown", - "id": "b6bff13f", + "id": "68999076", "metadata": {}, "source": [ "### Search and filter feature collection with pre-defined filters\n", @@ -155,7 +155,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7b9d374f", + "id": "4cb62faa", "metadata": {}, "outputs": [], "source": [ @@ -363,7 +363,7 @@ }, { "cell_type": "markdown", - "id": "45bd7b67", + "id": "41d6296c", "metadata": { "lines_to_next_cell": 0 }, @@ -382,7 +382,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d87cad2f", + "id": "f7e8ca10", "metadata": {}, "outputs": [], "source": [ @@ -418,7 +418,7 @@ }, { "cell_type": "markdown", - "id": "d8fabe5d", + "id": "bbe1f6e0", "metadata": { "lines_to_next_cell": 2 }, @@ -435,7 +435,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16182778", + "id": "b7097090", "metadata": { "lines_to_next_cell": 0 }, @@ -468,7 +468,7 @@ }, { "cell_type": "markdown", - "id": "cbb8ee58", + "id": "c50cec99", "metadata": {}, "source": [ "#### Find Laurussia topological plate and extract the reference features for investigation\n", @@ -484,7 +484,7 @@ { "cell_type": "code", "execution_count": null, - "id": "54860ec7", + "id": "614a6375", "metadata": {}, "outputs": [], "source": [ @@ -515,7 +515,7 @@ }, { "cell_type": "markdown", - "id": "b4fb459c", + "id": "dd52b0c4", "metadata": { "lines_to_next_cell": 0 }, @@ -530,7 +530,7 @@ { "cell_type": "code", "execution_count": null, - "id": "b49c55b4", + "id": "b752f2a9", "metadata": { "lines_to_next_cell": 0 }, @@ -561,7 +561,7 @@ }, { "cell_type": "markdown", - "id": "d5125a0d", + "id": "de192aa5", "metadata": {}, "source": [ "#### Topological reference map\n", @@ -577,7 +577,7 @@ { "cell_type": "code", "execution_count": null, - "id": "fd649426", + "id": "c2b8685f", "metadata": {}, "outputs": [], "source": [ @@ -595,7 +595,7 @@ }, { "cell_type": "markdown", - "id": "ac0440dc", + "id": "10a7e434", "metadata": { "lines_to_next_cell": 2 }, @@ -610,7 +610,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6835cda5", + "id": "b9ca5681", "metadata": {}, "outputs": [], "source": [ @@ -641,7 +641,7 @@ }, { "cell_type": "markdown", - "id": "eca479ff", + "id": "0be467ad", "metadata": {}, "source": [ "if you open the icosahedron_mesh_5.gpmlz file in GPlates, you will see the vertices of the icosahedron mesh,\n", @@ -655,7 +655,7 @@ }, { "cell_type": "markdown", - "id": "a512cac1", + "id": "4e0d6899", "metadata": { "lines_to_next_cell": 2 }, @@ -666,7 +666,7 @@ { "cell_type": "code", "execution_count": null, - "id": "22b6779a", + "id": "320a638d", "metadata": {}, "outputs": [], "source": [ @@ -694,7 +694,7 @@ }, { "cell_type": "markdown", - "id": "4937648e", + "id": "f75e060f", "metadata": {}, "source": [ "If you open icosahedron_vertices_in_region.gpmlz in GPlates, you will see the vertices of the icosahedron mesh that\n", @@ -709,7 +709,7 @@ }, { "cell_type": "markdown", - "id": "469147d9", + "id": "32031886", "metadata": { "lines_to_next_cell": 2 }, @@ -720,7 +720,7 @@ { "cell_type": "code", "execution_count": null, - "id": "194b9e5e", + "id": "e22c0d6f", "metadata": { "lines_to_next_cell": 0 }, @@ -766,7 +766,7 @@ }, { "cell_type": "markdown", - "id": "5f3e7dec", + "id": "f6dcfcc8", "metadata": {}, "source": [ "If you open icosahedron_vertices_within_australia.gpmlz in GPlates, you will see the vertices of the icosahedron mesh that\n", diff --git a/Notebooks/Examples/plot_map_with_pygmt.ipynb b/Notebooks/Examples/plot_map_with_pygmt.ipynb index bf1c073a..b4c17e18 100644 --- a/Notebooks/Examples/plot_map_with_pygmt.ipynb +++ b/Notebooks/Examples/plot_map_with_pygmt.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "a1c4908b", + "id": "32cfe43d", "metadata": {}, "source": [ "This notebook demonstrates how to use the PyGMT plotting maps in GPlately.\n", @@ -14,7 +14,7 @@ }, { "cell_type": "markdown", - "id": "fe001ac0", + "id": "4eb5307d", "metadata": {}, "source": [ "⚠️ This notebook is generated from plot_map_with_pygmt.py using the command `jupytext --to notebook Notebooks/Examples/plot_map_with_pygmt.py -o Notebooks/Examples/plot_map_with_pygmt.ipynb`.\n", @@ -25,7 +25,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26ee75f4", + "id": "52059e55", "metadata": {}, "outputs": [], "source": [ diff --git a/gplately/lib/reconstruct_by_topologies.py b/gplately/lib/reconstruct_by_topologies.py index f7c9c7ab..b032785d 100644 --- a/gplately/lib/reconstruct_by_topologies.py +++ b/gplately/lib/reconstruct_by_topologies.py @@ -14,7 +14,6 @@ # with this program; if not, write to Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. # -from abc import ABC, abstractmethod import logging import math import os @@ -29,401 +28,43 @@ logger = logging.getLogger("gplately") -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). +class ReconstructByTopologies(object): + """Reconstruct points using topologies and deactivates those that subduct at converging plate boundaries. + + This uses a new approach to collision detection of seed points + (to determine if they get subducted going forward in time and hence deactivated). + Previously we detected if the velocity delta of a seed point exceeded a threshold (when that point + transitioned from one plate to another) - and the point had to be close enough to the boundary. + Instead, we now take the plate containing a seed point at the current time and resolve it at the next time + (ie, current time plus time step) and see if its resolved shape still contains that seed point + (reconstructed to the next time). If it does then the seed point remains active (otherwise it's deactivated). + Note that the *same* plate is resolved at the current and next times. Usually when you resolve topologies + at a different time you will get different plates (eg, there could be plate merges or splits). + So, to resolve the *same* plate, some fenagling is required to ensure that a resolved plate boundary at the + current time can also be resolved at the next time. The end result is that any boundary around a plate that + converges will naturally deactivate seed points near that boundary (without having to use thresholds, etc). + And this doesn't require plate boundaries to be labelled as converging (ie, subduction zones). """ 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( # type: ignore - (curr_point,), - prev_location_velocity_stage_rotation, - 1, - pygplates.VelocityUnits.kms_per_my, - )[0] - curr_location_velocity = pygplates.calculate_velocities( # type: ignore - (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): - 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 - - -class _ReconstructByTopologies(ABC): - """Reconstruct geometries using topologies. Currently only points are supported.""" - - @abstractmethod - def _begin_reconstruction_impl(self): - """Called from 'begin_reconstruction()' and implemented in derived class. - - Derived class should call 'self._activate_deactivate_points_using_their_begin_end_times()' to introduce points - whose begin time is equal to (or older than) the initial time. - - Derived class can optionally do anything else needed before starting reconstruction over the time intervals. - """ - pass # This is an abstract method, no implementation here. - - @abstractmethod - def _reconstruct_to_next_time_impl(self): - """Called from 'reconstruct_to_next_time()' and implemented in derived class. - - Derived class should reconstruct the current active points in 'self.curr_points' to the next points 'self.next_points' from - the current time 'self.get_current_time()' to the next time 'self.get_current_time() + self.reconstruction_time_step'. - And collision detected for points should deactivate those points (ie, set None in 'self.next_points'). - - Derived class then should cycle the three point arrays in preparation for the next time interval. - Essentially do this: - self.prev_points, self.curr_points, self.next_points = self.curr_points, self.next_points, self.prev_points - - Derived class then should increment the current time index (to update the current time to the next time). - Essentially do this: - self.current_time_index += 1 - - Finally, derived class needs to call 'self._activate_deactivate_points_using_their_begin_end_times()' to introduce points - whose begin time is the new current time, and remove points whose end time is the new current time. - """ - pass # This is an abstract method, no implementation here. - - def __init__( - self, + plate_reconstruction, reconstruction_begin_time, reconstruction_end_time, reconstruction_time_interval, points, + *, point_begin_times: Union[np.ndarray, list, None] = None, point_end_times: Union[np.ndarray, list, None] = None, + continent_mask_filepath_format=None, ): + """ + continent_mask_filepath_format: str + The format of the file path of the continental mask grids that is converted to a + file path using ``continent_mask_filepath_format.format(time)``. + """ + + self.plate_reconstruction = plate_reconstruction # Set up an array of reconstruction times covering the reconstruction time span. self.reconstruction_begin_time = reconstruction_begin_time @@ -432,33 +73,33 @@ def __init__( 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 + self._reconstruction_time_step = -reconstruction_time_interval else: - self.reconstruction_time_step = reconstruction_time_interval + 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( + self._num_times = 1 + int( math.floor( 1e-6 + float(self.reconstruction_end_time - self.reconstruction_begin_time) - / self.reconstruction_time_step + / 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 # unless the reconstruction begin and end times are equal (in which case there'll be only one reconstruction snapshot). # 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 == 1 and self.reconstruction_end_time != self.reconstruction_begin_time ): - self.num_times = 2 - self.reconstruction_time_step = ( + 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._reconstruction_time_interval = math.fabs(self._reconstruction_time_step) - self.last_time_index = self.num_times - 1 + self._last_time_index = self._num_times - 1 self.points = points self.num_points = len(points) @@ -485,6 +126,8 @@ def __init__( "Length of 'point_end_times' must match length of 'points'." ) + self._continent_mask_filepath_format = continent_mask_filepath_format + def reconstruct(self): # Initialise the reconstruction. self.begin_reconstruction() @@ -495,562 +138,1161 @@ def reconstruct(self): while self.reconstruct_to_next_time(): pass - return self.get_active_current_points() + return self.get_current_active_points() def begin_reconstruction(self): - self.current_time_index = 0 + self._current_time_index = 0 + + # Store active and inactive points here. + # + # NOTE: Deactivated points will NOT get reset to None - we'll just ignore them. + self._all_current_points = [None] * self.num_points - # 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 + # A boolean array indicating which points (in 'self._all_current_points') are active. + self._point_is_currently_active = np.zeros(self.num_points, dtype=bool) # Each point can only get activated once (after deactivation it cannot be reactivated). - self.point_has_been_activated = np.zeros(self.num_points, dtype=bool) - self.num_activated_points = 0 + self._point_has_been_activated = np.zeros(self.num_points, dtype=bool) - # Call derived class implementation. - self._begin_reconstruction_impl() + # Activate any original points whose begin time is equal to (or older than) the newly updated current time. + # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. + self._activate_deactivate_points_using_their_begin_end_times() + + # Deactivate any newly actived points that are masked by the continents. + self._deactivate_continent_masked_points() def get_current_time_index(self): - return self.current_time_index + return self._current_time_index def get_current_time(self): return ( self.reconstruction_begin_time - + self.current_time_index * self.reconstruction_time_step + + self._current_time_index * self._reconstruction_time_step ) - def get_all_current_points(self): - return self.curr_points + def get_current_active_points(self): + """Return those points that are currently active.""" + return [ + self._all_current_points[point_index] + for point_index in self.get_current_active_point_indices() + ] - 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_current_active_point_indices(self): + """Return the indices of the currently active points (into the original points).""" + return np.where(self._point_is_currently_active)[0] 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: + 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): + if np.all(self._point_has_been_activated) and not np.any( + self._point_is_currently_active + ): return False - # Call derived class implementation. + # Call the main implementation. self._reconstruct_to_next_time_impl() - # We successfully reconstructed to the next time. - return True + # Move the current time to the next time. + self._current_time_index += 1 - def _activate_deactivate_points_using_their_begin_end_times(self): - current_time = self.get_current_time() + # Activate any original points whose begin time is equal to (or older than) the newly updated current time. + # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. + self._activate_deactivate_points_using_their_begin_end_times() - # 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 + # Deactivate any active points that are masked by the continents. + # + # This includes deactivating any newly activated points. + # + # Note: This is done at the newly updated current time. + self._deactivate_continent_masked_points() + # We successfully reconstructed to the next time. + return True -class _ReconstructByTopologiesImpl(_ReconstructByTopologies): - """Reconstruct geometries using topologies. Currently only points are supported.""" + def _reconstruct_to_next_time_impl(self): - use_plate_partitioner = False - """If the use_plate_partitioner is True then use pygplates.PlatePartitioner to partition points, - otherwise use faster points_in_polygons.find_polygons().""" + current_time = self.get_current_time() + # Positive/negative time step means reconstructing backward/forward in time. + next_time = current_time + self._reconstruction_time_step - def __init__( - self, - rotation_features_or_model, - topology_features, - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - *, - point_begin_times: Union[np.ndarray, list, None] = None, - point_end_times: Union[np.ndarray, list, None] = None, - point_plate_ids: Union[np.ndarray, list, None] = None, - detect_collisions=_DEFAULT_COLLISION, - ): - """ - rotation_features_or_model: Rotation model or feature collection(s), or list of features, or filename(s). + # Resolved the topologies at the current time. + curr_topological_snapshot = self.plate_reconstruction.topological_snapshot( + current_time, include_topological_slab_boundaries=False + ) - topology_features: Topology feature collection(s), or list of features, or filename(s) or any combination of those. + # Get the resolved networks and rigid plates. + curr_resolved_networks = curr_topological_snapshot.get_resolved_topologies( + pygplates.ResolveTopologyType.network # type: ignore + ) + curr_resolved_boundaries = curr_topological_snapshot.get_resolved_topologies( + pygplates.ResolveTopologyType.boundary # type: ignore + ) - detect_collisions: Collision detection function, or None. Defaults to _DEFAULT_COLLISION. - """ + # Get the currently active points and their indices (into the original points). + curr_active_point_indices = self.get_current_active_point_indices() + curr_active_points = [ + self._all_current_points[point_index] + for point_index in curr_active_point_indices + ] - super().__init__( - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - point_begin_times, - point_end_times, + # For each active current point find the resolved topology network or rigid plate containing it. + curr_active_points_spatial_tree = ( + _ptt.utils.points_spatial_tree.PointsSpatialTree(curr_active_points) + ) + curr_active_point_resolved_networks = ( + _ptt.utils.points_in_polygons.find_polygons_using_points_spatial_tree( + curr_active_points, + curr_active_points_spatial_tree, + [ + resolved_topology.get_resolved_boundary() + for resolved_topology in curr_resolved_networks + ], + curr_resolved_networks, + ) + ) + curr_active_point_resolved_boundaries = ( + _ptt.utils.points_in_polygons.find_polygons_using_points_spatial_tree( + curr_active_points, + curr_active_points_spatial_tree, + [ + resolved_topology.get_resolved_boundary() + for resolved_topology in curr_resolved_boundaries + ], + curr_resolved_boundaries, + ) ) - self.rotation_model = pygplates.RotationModel(rotation_features_or_model) + if len(curr_active_point_resolved_networks) != len( + curr_active_point_resolved_boundaries + ): + raise RuntimeError( + "Unexpected length mismatch of arrays of resolved topologies containing points." + ) - # Turn topology data into a list of features (if not already). - self.topology_features = pygplates.FeaturesFunctionArgument( - topology_features - ).get_features() + # Map of current resolved topologies to the active points contained within them (their point indices). + map_curr_resolved_topology_to_active_point_indices = {} - # 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. - if point_plate_ids is None: - self.point_plate_ids = np.zeros(self.num_points, dtype=int) - else: - # Make sure numpy array (if not already). - self.point_plate_ids = np.asarray(point_plate_ids) - if len(self.point_plate_ids) != self.num_points: - raise ValueError( - "Length of 'point_plate_ids' must match length of 'points'." + # Iterate over the resolved topologies containing all currently active points and get a list of the unique resolved topologies. + # + # Note: A resolved topology for a currently active point can be None if it fell outside all resolved topologies. + # In this case we'll just deactivate the point. Previously we would keep it active but just not reconstruct it + # (ie, the current and next positions would be the same). However the point might end up in weird locations. + # So it's best to just remove it. And this shouldn't happen very often at all (for topologies with global coverage). + for index in range(len(curr_active_points)): + + # Index into the active and inactive points 'self._all_current_points'. + curr_point_index = curr_active_point_indices[index] + + # If point is inside a resolved *network* then prefer that, otherwise see if it's inside a rigid plate. + # This is because deforming networks can overlay rigid plates. + curr_active_point_resolved_topology = curr_active_point_resolved_networks[ + index + ] + if curr_active_point_resolved_topology is None: + curr_active_point_resolved_topology = ( + curr_active_point_resolved_boundaries[index] ) - self.detect_collisions = detect_collisions - - def _begin_reconstruction_impl(self): - # Set up topology arrays (corresponding to active/inactive points at same indices as the point arrays). - 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 + # See if current active point fell outside all current resolved topologies. + if curr_active_point_resolved_topology is None: + # Current point is currently active but it fell outside all resolved topologies. + # So we deactivate it. + self._point_is_currently_active[curr_point_index] = False - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - # - # Note: This should be done before calling 'self._find_resolved_topologies_containing_points()' so that new - # points just appearing at their begin time will have associated topologies (that contain them). - self._activate_deactivate_points_using_their_begin_end_times() + # Continue to next active point. + continue - self._find_resolved_topologies_containing_points() + # If resolved topology has not been encountered yet then give it an empty list of point indices. + if ( + curr_active_point_resolved_topology + not in map_curr_resolved_topology_to_active_point_indices + ): + map_curr_resolved_topology_to_active_point_indices[ + curr_active_point_resolved_topology + ] = [] + + # Add the index of the currently active point to the resolved topology containing it. + map_curr_resolved_topology_to_active_point_indices[ + curr_active_point_resolved_topology + ].append(curr_point_index) + + # Prevent usage since might no longer be representative of the active points because + # we might've just deactivated some points that fell outside all resolved topologies. + del curr_active_points + del curr_active_point_indices + + # The current resolved topologies that contain active points. + curr_active_resolved_topologies = list( + map_curr_resolved_topology_to_active_point_indices.keys() + ) - def _reconstruct_to_next_time_impl(self): - # Cache stage rotations by plate ID. - # Use different dicts since using different rotation models and time steps, etc. - reconstruct_stage_rotation_dict = {} + # The *next* topology features and their topological section features. + next_topological_features = self._NextTopologicalFeatures(next_time) + # For each current resolved topology, create a corresponding *next* topological boundary feature and + # all the topological section features referenced by them. + for curr_resolved_topology in curr_active_resolved_topologies: + next_topological_features.add_current_topology(curr_resolved_topology) + + # Resolve the topologies to the *next* time step. + next_topological_snapshot = pygplates.TopologicalSnapshot( # type: ignore + next_topological_features.next_topological_section_features + + next_topological_features.next_topological_boundary_features, + self.plate_reconstruction.rotation_model, + next_time, + ) - current_time = self.get_current_time() + # First, map the *next* *boundary* features to their indices into 'next_topological_boundary_features'. + # + # This will also be the indices into 'curr_active_resolved_topologies' (due to order of insertion). + map_next_topology_boundary_feature_to_index = { + next_topology_feature: index + for index, next_topology_feature in enumerate( + next_topological_features.next_topological_boundary_features + ) + } + # + # Then, map each current resolved topology to its next resolved topology. + map_curr_to_next_resolved_topologies = {} + for ( + next_resolved_topology + ) in next_topological_snapshot.get_resolved_topologies(): + index = map_next_topology_boundary_feature_to_index.get( + next_resolved_topology.get_feature() + ) + # It's possible that the next topology could not get resolved at the next time step. + # This generally shouldn't happen though. + if index is not None: + curr_resolved_topology = curr_active_resolved_topologies[index] + map_curr_to_next_resolved_topologies[curr_resolved_topology] = ( + next_resolved_topology + ) - # 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 + # Iterate over all currently active points to reconstruct them to the next time step and + # see if they've been consumed by a convergent plate boundary. + # + # Note: Active points are grouped by the current resolved topology containing them. + for ( + curr_resolved_topology, + curr_resolved_topology_active_point_indices, + ) in map_curr_resolved_topology_to_active_point_indices.items(): + + # The next resolved topology associated with the current resolved topology. + next_resolved_topology = map_curr_to_next_resolved_topologies.get( + curr_resolved_topology + ) + if next_resolved_topology is None: + # The current topology could not be resolved to the next time step. + # So we deactivate all points that it contains. This generally shouldn't happen though. + self._point_is_currently_active[ + curr_resolved_topology_active_point_indices + ] = False 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, + # Get the boundary polygon of the next resolved topology. + # + # If the next resolved topology is consistent with the current resolved topology then we're fine, + # otherwise we need to replace the boundary of the next resolved topology with something more consistent. + next_resolved_topology_boundary = self._NextTopologicalBoundary( + curr_resolved_topology, + next_resolved_topology, + self.plate_reconstruction.rotation_model, + ).get_next_resolved_topology_boundary() + + # Iterate over the currently active points contained by the current resolved topology. + for point_index in curr_resolved_topology_active_point_indices: + + # Reconstruct the currently active point from its position at current time to its position at the next time step. + curr_active_point = self._all_current_points[point_index] + # This should not return None because the currently active point has previously been determined to be inside + # the current resolved topology (in the previous time step). + next_active_point = curr_resolved_topology.reconstruct_point( + curr_active_point, next_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 + # See if the location (of the current active point) reconstructed to the *next* time step + # is inside the current topology resolved to the *next* time step. + # + # If it is outside then it is subducting going forward in time (or consumed by a mid-ocean ridge going backward in time). + # It doesn't necessarily have to happen at a subduction zone (or mid-ocean ridge). It can be any part of a plate boundary + # that is convergent forward in time (or divergent forward in time; which is convergent backward in time). + if not next_resolved_topology_boundary.is_point_in_polygon( + next_active_point + ): + self._point_is_currently_active[point_index] = False + continue - # - # 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, + # The current active point was not consumed by a plate boundary. + # So update its position. + self._all_current_points[point_index] = next_active_point + + def _activate_deactivate_points_using_their_begin_end_times(self): + current_time = self.get_current_time() + + # Get indices of points that are not active and have never been activated. + curr_inactive_and_not_yet_activated_point_indices = np.where( + ~self._point_is_currently_active & ~self._point_has_been_activated + )[0] + # See which of those points we can activate. + activate_point_indices = curr_inactive_and_not_yet_activated_point_indices[ + ( + current_time + <= self.point_begin_times[ + curr_inactive_and_not_yet_activated_point_indices + ] + ) + & ( + current_time + >= self.point_end_times[ + curr_inactive_and_not_yet_activated_point_indices + ] + ) + ] + # Activate those points. + self._point_is_currently_active[activate_point_indices] = True + # Copy original point into currently active points array. + for point_index in activate_point_indices: + self._all_current_points[point_index] = self.points[point_index] + # Mark those points as having been activated. + self._point_has_been_activated[activate_point_indices] = True + + # Get indices of points that are active (and hence must have been previously activated). + curr_active_point_indices = np.where(self._point_is_currently_active)[0] + # See which of those points we can deactivate. + deactivate_point_indices = curr_active_point_indices[ + ~( + (current_time <= self.point_begin_times[curr_active_point_indices]) + & (current_time >= self.point_end_times[curr_active_point_indices]) + ) + ] + # Deactivate those points. + self._point_is_currently_active[deactivate_point_indices] = False + # Leave deactivated points in the currently active points array (we'll just ignore them). + + def _deactivate_continent_masked_points(self): + # If no continental collision detection requested then just return. + if self._continent_mask_filepath_format is None: + return + + # Get the currently active points and their indices (into the original points). + curr_active_point_indices = self.get_current_active_point_indices() + if len(curr_active_point_indices) == 0: + # There are no active points so nothing to do. + return + curr_active_points = [ + self._all_current_points[point_index] + for point_index in curr_active_point_indices + ] + + # Convert points to lat/lon. + curr_points_lat = np.empty(len(curr_active_points)) + curr_points_lon = np.empty(len(curr_active_points)) + for index, point in enumerate(curr_active_points): + point_lat, point_lon = point.to_lat_lon() + curr_points_lat[index] = point_lat + curr_points_lon[index] = point_lon + + # Read the continent mask at the current time. + continent_mask_filepath = self._continent_mask_filepath_format.format( + self.get_current_time() ) - self.prev_resolved_plate_boundaries, self.curr_resolved_plate_boundaries = ( - self.curr_resolved_plate_boundaries, - self.prev_resolved_plate_boundaries, + gridZ, gridX, gridY = read_netcdf_grid( + continent_mask_filepath, return_grids=True ) - # - # Move the current time to the next time. - self.current_time_index += 1 - current_time = self.get_current_time() + ni, nj = gridZ.shape + xmin = np.nanmin(gridX) + xmax = np.nanmax(gridX) + ymin = np.nanmin(gridY) + ymax = np.nanmax(gridY) + + # Sample continent mask grid, which is one over continents and zero over oceans. + points_i = (ni - 1) * ((curr_points_lat - ymin) / (ymax - ymin)) + points_j = (nj - 1) * ((curr_points_lon - xmin) / (xmax - xmin)) + points_i_uint = np.rint(points_i).astype(np.uint) + points_j_uint = np.rint(points_j).astype(np.uint) + try: + mask_values = gridZ[points_i_uint, points_j_uint] + except IndexError: + points_i = np.clip(np.rint(points_i), 0, ni - 1).astype(np.int_) + points_j = np.clip(np.rint(points_j), 0, nj - 1).astype(np.int_) + mask_values = gridZ[points_i, points_j] - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - # - # Note: This should be done before calling 'self._find_resolved_topologies_containing_points()' so that new - # points just appearing at their begin time will have associated topologies (that contain them). - self._activate_deactivate_points_using_their_begin_end_times() + # Deactivate any points that sampled inside the continent mask. + self._point_is_currently_active[ + curr_active_point_indices[np.where(mask_values >= 0.5)[0]] + ] = False - 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 + class _NextTopologicalFeatures(object): + """The *next* topology boundary features and their topological section features. - # 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] + These are essentially the same as the *current* topology features (and their topological section features) but + with their valid time periods modified (if needed) to include the *next* time. + """ - # 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], + def __init__(self, next_time): + + self.next_time = next_time + # All topology *boundary* features and their referenced topological *section* features ready to be resolved at the next time. + # + # Note: Topological *section* features can include topological lines. + self.next_topological_boundary_features = [] + self.next_topological_section_features = [] + + # Map of each current topological section feature to its corresponding *next* topological section feature. + self._map_curr_to_next_topological_section_feature = {} + + def add_current_topology(self, curr_resolved_topology): + + # Iterate over the boundary sub-segments of the current resolved topology and create a new topological + # section feature for each one (if its topological section feature hasn't been encountered yet). + next_boundary_section_property_values = [] + for ( + curr_boundary_sub_segment + ) in curr_resolved_topology.get_boundary_sub_segments(): # type: ignore + + curr_boundary_section_feature = curr_boundary_sub_segment.get_feature() + # If the current boundary section feature hasn't been encountered yet then create an associated + # *next* boundary section feature that is either the same as the current boundary section feature + # (if it's still valid at the *next* time) or a clone of the current boundary section feature with + # the valid time range modified to be valid at the *next* time. + if ( + curr_boundary_section_feature + not in self._map_curr_to_next_topological_section_feature ): - # 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 + curr_boundary_section_recon_geom = ( + curr_boundary_sub_segment.get_topological_section() + ) + # Boundary section is either a ReconstructedFeatureGeometry or a ResolvedTopologicalLine. + # + # If it's a ResolvedTopologicalLine then we need to iterate through its sub-segments and create a + # new line section for each one, and finally create a new topological line feature from them. + if isinstance( + curr_boundary_section_recon_geom, pygplates.ResolvedTopologicalLine # type: ignore + ): - def _find_resolved_topologies_containing_points(self): - current_time = self.get_current_time() + # Iterate over the sub-segments of the current resolved topological line and create a new line + # section feature for each one (if its topological section feature hasn't been encountered yet). + next_line_section_property_values = [] + for ( + curr_line_sub_segment + ) in curr_boundary_section_recon_geom.get_line_sub_segments(): # type: ignore + curr_line_section_feature = ( + curr_line_sub_segment.get_feature() + ) + # If the current line section feature hasn't been encountered yet then create an associated + # *next* line section feature that is either the same as the current line section feature + # (if it's still valid at the *next* time) or a clone of the current line section feature with + # the valid time range modified to be valid at the *next* time. + if ( + curr_line_section_feature + not in self._map_curr_to_next_topological_section_feature + ): + # Line section must be a ReconstructedFeatureGeometry (ie, can't be ResolvedTopologicalLine) + # since a topological line cannot have sections that are also topological lines. + + # If feature is not valid at 'next_time' then clone the feature and set a valid time that includes 'next_time'. + if not curr_line_section_feature.is_valid_at_time( + self.next_time + ): + next_line_section_feature = ( + curr_line_section_feature.clone() + ) + next_line_section_feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + else: + next_line_section_feature = ( + curr_line_section_feature + ) + + self.next_topological_section_features.append( + next_line_section_feature + ) + + # Map the current to the next (line section feature). + self._map_curr_to_next_topological_section_feature[ + curr_line_section_feature + ] = next_line_section_feature + + # Get the next line section feature associated with the current one. + next_line_section_feature = ( + self._map_curr_to_next_topological_section_feature[ + curr_line_section_feature + ] + ) + + # The property name of the geometry referenced by the next line section. + next_line_section_geometry_property_name = ( + curr_line_sub_segment.get_topological_section() + .get_property() + .get_name() + ) + # Whether to reverse the geometry referenced by the next line section. + next_line_section_reverse_order = ( + curr_line_sub_segment.was_geometry_reversed_in_topology() + ) + + # Create the *next* line section associated with the current one. + next_line_section_property_value = pygplates.GpmlTopologicalSection.create( # type: ignore + next_line_section_feature, + geometry_property_name=next_line_section_geometry_property_name, + reverse_order=next_line_section_reverse_order, + topological_geometry_type=pygplates.GpmlTopologicalLine, # type: ignore + ) + if not next_line_section_property_value: + raise RuntimeError( + "Unable to create topological section for a topological line." + ) + next_line_section_property_values.append( + next_line_section_property_value + ) + + # Create a topological line feature that enables the current topological line + # to be resolved to the *next* time step. + # + # Note: 'valid_time' argument is not specified which means valid for all time. + # This works for us since we only need it to be valid at the *next* time step. + next_line_feature = pygplates.Feature( # type: ignore + pygplates.FeatureType.gpml_unclassified_feature # type: ignore + ) + next_line_feature.set_topological_geometry( + pygplates.GpmlTopologicalLine( # type: ignore + next_line_section_property_values + ), + # Use the same geometry property name as the current boundary sub-segment so + # that it can be found by the topological polygon(s) that will reference it... + curr_boundary_sub_segment.get_topological_section() + .get_property() + .get_name(), + ) + + next_boundary_section_feature = next_line_feature + + else: # boundary section is a ReconstructedFeatureGeometry ... + + # If feature is not valid at 'next_time' then clone the feature and set a valid time that includes 'next_time'. + if not curr_boundary_section_feature.is_valid_at_time( + self.next_time + ): + next_boundary_section_feature = ( + curr_boundary_section_feature.clone() + ) + next_boundary_section_feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + else: + next_boundary_section_feature = ( + curr_boundary_section_feature + ) + + self.next_topological_section_features.append( + next_boundary_section_feature + ) - # Resolve the plate polygons for the current time. - resolved_topologies = [] - pygplates.resolve_topologies( # type: ignore - self.topology_features, - self.rotation_model, - resolved_topologies, - current_time, - ) + # Map the current to the next (boundary section feature). + self._map_curr_to_next_topological_section_feature[ + curr_boundary_section_feature + ] = next_boundary_section_feature - if _ReconstructByTopologiesImpl.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) # type: ignore - 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, + # Get the next boundary section feature associated with the current one. + next_boundary_section_feature = ( + self._map_curr_to_next_topological_section_feature[ + curr_boundary_section_feature + ] ) - ) - # 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 + # The property name of the geometry referenced by the next boundary section. + next_boundary_section_geometry_property_name = ( + curr_boundary_sub_segment.get_topological_section() + .get_property() + .get_name() + ) + # Whether to reverse the geometry referenced by the next boundary section. + next_boundary_section_reverse_order = ( + curr_boundary_sub_segment.was_geometry_reversed_in_topology() + ) - # Find the plate id of the polygon that contains 'curr_point'. - if _ReconstructByTopologiesImpl.use_plate_partitioner: - curr_polygon = plate_partitioner.partition_point(curr_point) # type: ignore - else: - curr_polygon = resolved_topologies_containing_curr_valid_points[ # type: ignore - # Index back into 'curr_valid_points' and hence also into - # 'resolved_topologies_containing_curr_valid_points'. - curr_valid_points_indices[point_index] # type: ignore - ] # type: ignore - 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 + # Create the *next* boundary section associated with the current one. + next_boundary_section_property_value = pygplates.GpmlTopologicalSection.create( # type: ignore + next_boundary_section_feature, + geometry_property_name=next_boundary_section_geometry_property_name, + reverse_order=next_boundary_section_reverse_order, + topological_geometry_type=pygplates.GpmlTopologicalPolygon, # type: ignore + ) + if not next_boundary_section_property_value: + raise RuntimeError( + "Unable to create topological section for a topological polygon." + ) + next_boundary_section_property_values.append( + next_boundary_section_property_value + ) - # 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() + # Create a topological boundary feature that enables the current topological boundary + # to be resolved to the *next* time step. + # + # Note: If the current topology is a deforming network we still create a topological closed plate boundary + # (which is normally used for rigid plates) because we only need to detect if a seed point + # reconstructed to the next time step is contained within the network's *boundary*. + # + # Note: 'valid_time' argument is not specified which means valid for all time. + # This works for us since we only need it to be valid at the *next* time step. + next_boundary_feature = pygplates.Feature.create_topological_feature( # type: ignore + pygplates.FeatureType.gpml_topological_closed_plate_boundary, # type: ignore + pygplates.GpmlTopologicalPolygon( # type: ignore + next_boundary_section_property_values + ), ) + self.next_topological_boundary_features.append(next_boundary_feature) + + class _NextTopologicalBoundary(object): + """The next resolved topology boundary. + + If the next resolved topology is NOT consistent with the current resolved topology then we need to replace + replace any inconsistent boundary sub-segments of the next resolved topology with the equivalent ones in the + current resolved topology (but moved to the next time). For each inconsistent boundary sub-segment, this involves + taking the boundary sub-segment of the *current* resolved topology and moving it to the *next* time, and using that + (instead of the boundary sub-segment of the *next* resolved topology). This is not as ideal as actually resolving the + all the *unclipped* boundary sections of the current topology at the *next* time, which would produce a more accurate + plate boundary for the next resolved topology. But the inconsistency forces us to take an alternative approach. + However this case shouldn't happen very often. + + The next resolved topology is not consistent with the current resolved topology if its boundary sections don't have + the same number of intersections with neighbouring sections. For example, a boundary section of the current resolved topology + intersects its neighbour once (as expected) but the same boundary section of the next resolved topology (which is really + just the same topology resolved to the next time) intersects the same neighbour twice. In this case, the next resolved topology + would likely have an unexpectedly different plate boundary shape than the current resolved topology, which might deactivate + a bunch of points that it shouldn't. This happens because pyGPlates only considers the first intersection (if there are two or more). + Usually the topological model is built such that each neighbouring section only intersects once. However, when you resolve + that same topology at a different time (even if it's just 1 Myr different) then the intersections may not be what the + model builder intended (especially if the current topology's time period does not include the *next* time). + """ + # Distance from a polyline-polyline intersection such that if we clipped that much distance off the polylines + # around that intersection then the clipped polylines are unlikely to intersect (unless they're almost parallel). + # + # We are a little conversative on this (ie, a little big) to ensure we don't get too close. + INTERSECTION_THRESHOLD_RADIANS = 1e-4 # approx 600 metres -class _ReconstructByTopologicalModelImpl(_ReconstructByTopologies): - """Similar to ``_ReconstructByTopologiesImpl`` except uses the ``pygplates.TopologicalModel`` class to reconstruct seed points. + # Maximum number of intersections to test between two intersecting polylines. + # + # If two polylines reach the max number of intersections then we've likely detected a partial overlap + # between the two polylines (where a segment from each polyline are on top of each other and + # hence essentially have an infinite number of intersections). + MAX_NUM_INTERSECTIONS = 5 - This is currently just a transition towards using ``pygplates.TopologicalModel``. - But note that this still uses Python code, like in class ``_DefaultCollision``, to do the collision detection - (as opposed to collision detection built into pyGPlates, which needs to be updated to better support GPlately). - Also it's not as fast as ``_ReconstructByTopologiesImpl`` which has faster point-in-polygon testing (since it uses a spatial tree). + def __init__( + self, curr_resolved_topology, next_resolved_topology, rotation_model + ): - So, for the time being it's probably still better to use class ``_ReconstructByTopologiesImpl`` instead. - """ + self.curr_resolved_topology = curr_resolved_topology + self.next_resolved_topology = next_resolved_topology - class CollisionDelegator(pygplates.ReconstructedGeometryTimeSpan.DeactivatePoints): - """Delegates collision detection from pyGPlates (pygplates.TopologicalModel) to the collision classes accessed by ``_ReconstructByTopologies`` - (like class ``_DefaultCollision``). + self.curr_boundary_sub_segments = ( + curr_resolved_topology.get_boundary_sub_segments() + ) + self.next_boundary_sub_segments = ( + next_resolved_topology.get_boundary_sub_segments() + ) - Ultimately there should be less (or no) need to delegate once the default pyGPlates collision handler (``DefaultDeactivatePoints``) supports the - funcionality of classes like ``_DefaultCollision``. Although ``_ContinentCollision`` might still need to be handled outside pyGPlates (ie, delegated to). - """ + if len(self.curr_boundary_sub_segments) != len( + self.next_boundary_sub_segments + ): + raise RuntimeError( + "Expected current and next topologies have the same number of boundary sub-segments." + ) + self.num_boundary_sub_segments = len(self.curr_boundary_sub_segments) - def __init__( - self, detect_collisions, rotation_model, reconstruction_time_interval - ): - super().__init__() - self.detect_collisions = detect_collisions + self.curr_time = curr_resolved_topology.get_reconstruction_time() + self.next_time = next_resolved_topology.get_reconstruction_time() self.rotation_model = rotation_model - self.reconstruction_time_interval = reconstruction_time_interval - - def deactivate( - self, - prev_point, - prev_location, - prev_time, - curr_point, - curr_location, - curr_time, + + def get_next_resolved_topology_boundary(self): + + # If the current and next resolved topologies are consistent then + # just return the resolved boundary of the next resolved topology. + inconsistent_boundary_sub_segment_indices = ( + self._get_inconsistent_boundary_sub_segments() + ) + if not inconsistent_boundary_sub_segment_indices: + return self.next_resolved_topology.get_resolved_boundary() + + # Instead of returning the boundary of the next resolved topology, replace any inconsistent boundary sub-segments + # of the next resolved topology with the equivalent ones in the current resolved topology (but moved to the next time). + return self._get_next_boundary_from_consistent_boundary_sub_segments( + inconsistent_boundary_sub_segment_indices + ) + + def _get_next_boundary_from_consistent_boundary_sub_segments( + self, inconsistent_boundary_sub_segment_indices ): - # Get plate IDs of resolved topology containing previous and current point. + + # + # Instead of the resolved boundary of the *next* resolved topology we'll replace the inconsistent boundary sub-segments + # of the *next* resolved topology with the equivalent ones in the *current* resolved topology (but moved to the next time). # - # Note that could be None, so the collision detection needs to handle that. - prev_plate_id = None - curr_plate_id = None - - # See if previous point is located in a resolved rigid plate or deforming network (or neither). - prev_resolved_topology = prev_location.located_in_resolved_boundary() - if not prev_resolved_topology: - prev_resolved_topology = prev_location.located_in_resolved_network() - if prev_resolved_topology: - prev_plate_id = ( - prev_resolved_topology.get_feature().get_reconstruction_plate_id() + + # Iterate over the next boundary sub-segments and join them to form a polygon boundary. + next_boundary_polygon_points = [] + for boundary_sub_segment_index in range(self.num_boundary_sub_segments): + + # If the boundary sub-segment is consistent (ie, not inconsistent) then just add the + # (potentially reversed) points of the boundary sub-segment of the *next* resolved topology. + if ( + boundary_sub_segment_index + not in inconsistent_boundary_sub_segment_indices + ): + next_boundary_sub_segment = self.next_boundary_sub_segments[ + boundary_sub_segment_index + ] + next_boundary_polygon_points.extend( + reversed( + next_boundary_sub_segment.get_resolved_geometry_points() + ) + if next_boundary_sub_segment.was_geometry_reversed_in_topology() + else next_boundary_sub_segment.get_resolved_geometry_points() + ) + + continue + + # + # The boundary sub-segment is inconsistent. + # + # So take the boundary sub-segment of the *current* resolved topology and move it to the *next* time, + # and use that (instead of the boundary sub-segment of the *next* resolved topology). + # + + curr_boundary_sub_segment = self.curr_boundary_sub_segments[ + boundary_sub_segment_index + ] + curr_boundary_sub_segment_reversal = ( + curr_boundary_sub_segment.was_geometry_reversed_in_topology() ) - # See if current point is located in a resolved rigid plate or deforming network (or neither). - curr_resolved_topology = curr_location.located_in_resolved_boundary() - if not curr_resolved_topology: - curr_resolved_topology = curr_location.located_in_resolved_network() - if curr_resolved_topology: - curr_plate_id = ( - curr_resolved_topology.get_feature().get_reconstruction_plate_id() + # See if the boundary sub-segment is from a topological *line*. + curr_boundary_sub_sub_segments = ( + curr_boundary_sub_segment.get_sub_segments() ) + if not curr_boundary_sub_sub_segments: + # + # Boundary sub-segment is NOT from a topological *line*. + # + curr_boundary_sub_segment_feature = ( + curr_boundary_sub_segment.get_resolved_feature() + ) - # Delegate to the Python implementation of collision detection (eg, classes '_DefaultCollision' and '_ContinentCollision'). - return self.detect_collisions( - self.rotation_model, - curr_time, - self.reconstruction_time_interval, - prev_point, - curr_point, - prev_plate_id, - prev_resolved_topology, - curr_plate_id, - curr_resolved_topology, - ) + # The sub-segment resolved feature is already reconstructed to 'curr_time'. + # So we need to reverse reconstruct it to present day (before we can reconstruct it to 'next_time'). + # + # Note: We don't need to clone this feature before modifying it because it was already essentially cloned + # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. + pygplates.reverse_reconstruct( # type: ignore + curr_boundary_sub_segment_feature, + self.rotation_model, + self.curr_time, + ) - def __init__( - self, - rotation_features_or_model, - topology_features, - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - *, - point_begin_times: Union[np.ndarray, list, None] = None, - point_end_times: Union[np.ndarray, list, None] = None, - detect_collisions=_DEFAULT_COLLISION, - ): - """ - rotation_features_or_model: Rotation model or feature collection(s), or list of features, or filename(s). + # The resolved sub-segment feature is valid at the 'curr_time' but not necessarily at 'next_time'. + # If not valid at 'next_time' then set a valid time that includes 'next_time'. + if not curr_boundary_sub_segment_feature.is_valid_at_time( + self.next_time + ): + curr_boundary_sub_segment_feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + + # Reconstruct to 'next_time' (from present day). + next_boundary_sub_segment_reconstructed_geometries = ( + pygplates.ReconstructSnapshot( # type: ignore + curr_boundary_sub_segment_feature, + self.rotation_model, + self.next_time, + ).get_reconstructed_geometries() + ) - topology_features: Topology feature collection(s), or list of features, or filename(s) or any combination of those. + # Should only be one reconstructed geometry. + if len(next_boundary_sub_segment_reconstructed_geometries) != 1: + raise RuntimeError("Expected a single reconstructed geometry.") + next_boundary_sub_segment_reconstructed_geometry = ( + next_boundary_sub_segment_reconstructed_geometries[0] + ) - detect_collisions: Collision detection function, or None. Defaults to _DEFAULT_COLLISION. - """ + # Add the (potentially reversed) points of the sub-segment to the next polygon boundary. + next_boundary_polygon_points.extend( + reversed( + next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) + if curr_boundary_sub_segment_reversal + else next_boundary_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) - super().__init__( - reconstruction_begin_time, - reconstruction_end_time, - reconstruction_time_interval, - points, - point_begin_times, - point_end_times, - ) + continue - self.topological_model = pygplates.TopologicalModel( - topology_features, - rotation_features_or_model, - # Only really need to cache 2 topological snapshots since we progressively move through time (and hence never return to previous times). - # Might only need to cache 1 topological snapshot (not sure). Best to be safe with 2... - topological_snapshot_cache_size=2, - ) + # + # The boundary sub-segment IS from a topological *line*, so we need to iterate over its sub-segments. + # This is because we need features that are *reconstructable*. + # Topological lines are not reconstructable (because they're topologies linking reconstructable features). + # Only the sub-segments of topological lines are reconstructable. + # + # Each sub-segment of the resolved topological *line* will contribute to the boundary sub-segment. + # + # Note: Only a sub-section of the full resolved topological *line* will contribute + # to the boundary of the current resolved topological boundary/network. + # Hence not all the sub-segments of the resolved topological *line* will contribute + # to the boundary of the current resolved topological boundary/network. + # - self.detect_collisions = self.CollisionDelegator( - detect_collisions, - self.topological_model.get_rotation_model(), - self.reconstruction_time_interval, - ) + # Collect the sub-sub-segment features to use for the next boundary. + # + # Note: We need to reverse the order of sub-sub-segments if the resolved topological *line* is reversed + # in the current resolved boundary. + curr_boundary_sub_sub_segment_features = [ + sub_sub_segment.get_resolved_feature() + for sub_sub_segment in ( + reversed(curr_boundary_sub_sub_segments) + if curr_boundary_sub_segment_reversal + else curr_boundary_sub_sub_segments + ) + ] + # If a sub-sub-segment is reversed in the resolved topological *line* which is, in turn, reversed in the + # current resolved boundary then the sub-sub-segment is NOT reversed in the current resolved boundary. + curr_boundary_sub_sub_segment_reversals = [ + curr_boundary_sub_segment_reversal + ^ sub_sub_segment.was_geometry_reversed_in_topology() + for sub_sub_segment in ( + reversed(curr_boundary_sub_sub_segments) + if curr_boundary_sub_segment_reversal + else curr_boundary_sub_sub_segments + ) + ] - def _begin_reconstruction_impl(self): - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - self._activate_deactivate_points_using_their_begin_end_times() + # The sub-segment resolved features are already reconstructed to 'curr_time'. + # So we need to reverse reconstruct them to present day (before we can reconstruct them to 'next_time'). + # + # Note: We don't need to clone these features before modifying them because they were already essentially cloned + # when 'ResolvedTopologicalSubSegment.get_resolved_feature()' was called. + pygplates.reverse_reconstruct( # type: ignore + curr_boundary_sub_sub_segment_features, + self.rotation_model, + self.curr_time, + ) - def _reconstruct_to_next_time_impl(self): + # The resolved sub-segment features are valid at the 'curr_time' but not necessarily at 'next_time'. + # If not valid at 'next_time' then set a valid time that includes 'next_time'. + for feature in curr_boundary_sub_sub_segment_features: + if not feature.is_valid_at_time(self.next_time): + feature.set_valid_time( + pygplates.GeoTimeInstant.create_distant_past(), 0.0 # type: ignore + ) + + # Reconstruct to 'next_time' (from present day). + next_boundary_sub_sub_segment_reconstructed_geometries = pygplates.ReconstructSnapshot( # type: ignore + curr_boundary_sub_sub_segment_features, + self.rotation_model, + self.next_time, + ).get_reconstructed_geometries( + # Make sure the order is retained... + same_order_as_reconstructable_features=True + ) - current_time = self.get_current_time() - next_time = current_time + self.reconstruction_time_step - - # Find the active points. - # These will be reconstructed to the next time step. - current_active_points = [] - current_active_point_indices = [] - for point_index in range(self.num_points): - current_point = self.curr_points[point_index] - if current_point is None: - # Current point is not currently active. - # So we cannot reconstruct to next time. - self.next_points[point_index] = None - continue + # The order and length of reconstructed geometries should match that of the features. + if len(next_boundary_sub_sub_segment_reconstructed_geometries) != len( + curr_boundary_sub_sub_segment_features + ): + raise RuntimeError( + "Expected number of reconstructed geometries to match number of features." + ) - # Keep track of the currently active points. - # And their original indices into ALL points (active and inactive). - current_active_points.append(current_point) - current_active_point_indices.append(point_index) - - # Reconstruct active points to the next time step. - reconstructed_time_span = self.topological_model.reconstruct_geometry( - current_active_points, - initial_time=current_time, - youngest_time=next_time, - time_increment=self.reconstruction_time_interval, # must be positive - deactivate_points=self.detect_collisions, - ) + # Add the (potentially reversed) points of the sub-segments to the next polygon boundary. + for ( + feature_index, + sub_sub_segment_reconstructed_geometry, + ) in enumerate(next_boundary_sub_sub_segment_reconstructed_geometries): + next_boundary_polygon_points.extend( + reversed( + sub_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) + if curr_boundary_sub_sub_segment_reversals[feature_index] + else sub_sub_segment_reconstructed_geometry.get_reconstructed_geometry_points() + ) - # Store the next points back to their original locations in ALL points (active and inactive). - next_points = reconstructed_time_span.get_geometry_points( - next_time, return_inactive_points=True - ) - for next_point_index, next_point in enumerate(next_points): - # Get index into ALL points (active and inactive). - point_index = current_active_point_indices[next_point_index] - # Update next point (note that this can be None if a collision was detected). - self.next_points[point_index] = next_point + # Join the (potentially reversed) points of the boundary sub-segments and join them together into a polygon boundary. + # Each sub-segment is either from the *current* or *next* resolved topology. + return pygplates.PolygonOnSphere( # type: ignore + next_boundary_polygon_points + ) - # - # 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, - ) - # - # Move the current time to the next time. - self.current_time_index += 1 + def _get_inconsistent_boundary_sub_segments(self): + """Returns a set of indices of any boundary sub-segments that are inconsistent, otherwise None.""" - # Activate any original points whose begin time is equal to (or older than) the newly updated current time. - # Deactivate any activated points whose end time is equal to (or younger than) the newly updated current time. - self._activate_deactivate_points_using_their_begin_end_times() + # If there's only 1 sub-segment then it cannot intersect with itself (but it can still be converted to a polygon). + # + # Note: If there's 2 sub-segments then they can intersect each other twice (to form a closed polygon). + # We will still process them to make sure they stay consistent (same number of intersections at current and next times). + if self.num_boundary_sub_segments <= 1: + return None + + inconsistent_boundary_sub_segment_indices = None + + # Compare number of neighbour intersections for the current and next boundary sub-segments. + for boundary_sub_segment_index in range(self.num_boundary_sub_segments): + curr_boundary_sub_segment = self.curr_boundary_sub_segments[ + boundary_sub_segment_index + ] + next_boundary_sub_segment = self.next_boundary_sub_segments[ + boundary_sub_segment_index + ] + + prev_boundary_sub_segment_index = boundary_sub_segment_index - 1 + if prev_boundary_sub_segment_index < 0: + prev_boundary_sub_segment_index += self.num_boundary_sub_segments + + # Previous neighbour sub-segment (for the current and next resolved topologies). + prev_curr_boundary_sub_segment = self.curr_boundary_sub_segments[ + prev_boundary_sub_segment_index + ] + prev_next_boundary_sub_segment = self.next_boundary_sub_segments[ + prev_boundary_sub_segment_index + ] + + # Number of neighbour intersections of sub-segment (for the current and next resolved topologies). + curr_num_intersections = self._get_num_intersections( + prev_curr_boundary_sub_segment.get_topological_section_geometry(), + curr_boundary_sub_segment.get_topological_section_geometry(), + ) + next_num_intersections = self._get_num_intersections( + prev_next_boundary_sub_segment.get_topological_section_geometry(), + next_boundary_sub_segment.get_topological_section_geometry(), + ) + + # If the number of intersections differ then the current and next resolved topologies are inconsistent. + if curr_num_intersections != next_num_intersections: + # We only create a set when we actually need one + # (most often we won't since most resolved topologies are consistent). + if inconsistent_boundary_sub_segment_indices is None: + inconsistent_boundary_sub_segment_indices = set() + # Add the previous and current sub-segments affected by the intersection(s). + inconsistent_boundary_sub_segment_indices.add( + prev_boundary_sub_segment_index + ) + inconsistent_boundary_sub_segment_indices.add( + boundary_sub_segment_index + ) + + return inconsistent_boundary_sub_segment_indices + + def _get_num_intersections(self, polyline1, polyline2, num_intersections=0): + + distance_tuple = pygplates.GeometryOnSphere.distance( # type: ignore + polyline1, + polyline2, + # We're only detecting zero distance (an intersection), so as an optimisation we + # don't need to calculate an accurate non-zero distance (instead just return None)... + distance_threshold_radians=self.INTERSECTION_THRESHOLD_RADIANS, + return_closest_positions=True, + return_closest_indices=True, + ) + + if distance_tuple is not None: + ( + distance, + polyline1_intersection, + polyline2_intersection, + polyline1_segment_index, + polyline2_segment_index, + ) = distance_tuple + + if distance == 0.0: + # We found an intersection (between 'polyline1' and 'polyline2'). + num_intersections += 1 + + # If we reached the max number of intersections then we've likely detected a partial overlap + # between the two polylines (where a segment from each polyline are on top of each other and + # hence essentially have an infinite number of intersections). + # + # In this case we'll just return the maximum number of intersections. This means if both the + # polylines overlap at both the current and next times then they'll return the same number + # of intersections (max) and compare equal and hence be considered consistent. + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + + # Split 'polyline1' into two polylines at its intersection point. + polyline1a = self._pre_split_polyline( + polyline1, polyline1_segment_index, polyline1_intersection + ) + polyline1b = self._post_split_polyline( + polyline1, polyline1_segment_index, polyline1_intersection + ) + + # Split 'polyline2' into two polylines at its intersection point. + polyline2a = self._pre_split_polyline( + polyline2, polyline2_segment_index, polyline2_intersection + ) + polyline2b = self._post_split_polyline( + polyline2, polyline2_segment_index, polyline2_intersection + ) + + # Test for intersections among the four combinations of split polylines. + # + # Note: It's possible there may be less than four split polylines if + # either (or both) polylines were intersected at one of their end points. + if polyline1a: + if polyline2a: + num_intersections = self._get_num_intersections( + polyline1a, polyline2a, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + if polyline2b: + num_intersections = self._get_num_intersections( + polyline1a, polyline2b, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + if polyline1b: + if polyline2a: + num_intersections = self._get_num_intersections( + polyline1b, polyline2a, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + if polyline2b: + num_intersections = self._get_num_intersections( + polyline1b, polyline2b, num_intersections + ) + if num_intersections == self.MAX_NUM_INTERSECTIONS: + return num_intersections + + return num_intersections + + def _pre_split_polyline( + self, polyline, polyline_segment_index, polyline_intersection + ): + # Get the points *before* the intersection (including the intersection). + # This will have at least two points (needed for creating a polyline). + pre_split_points = list(polyline[: polyline_segment_index + 1]) + pre_split_points.append(polyline_intersection) + + # Remove any zero length segments next to the intersected end of the polyline. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + while True: + pre_split_last_segment = pygplates.GreatCircleArc( # type: ignore + pre_split_points[-2], + pre_split_points[-1], + ) + if not pre_split_last_segment.is_zero_length(): + break + + # Remove the last segment (by removing last point). + del pre_split_points[-1] + if len(pre_split_points) < 2: + # All segments before the intersection are zero length. + # So there's no pre-split polyline. + return None + + # We've encountered a non-zero length segment. Let's shorten it such that it will no longer intersect the intersection point. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + if ( + self.INTERSECTION_THRESHOLD_RADIANS + < pre_split_last_segment.get_arc_length() + ): + # Last segment is *longer* than the intersection threshold, so shorten it by that amount. + # Replace last point with rotated point. + pre_split_points[-1] = ( + pygplates.FiniteRotation( # type: ignore + pre_split_last_segment.get_rotation_axis(), + -self.INTERSECTION_THRESHOLD_RADIANS, # negative rotates away from segment *end* point + ) + * pre_split_points[-1] + ) + else: + # Last segment is *shorter* than the intersection threshold. But we know it's not a zero length segment, + # so it should be far enough away from the intersection point. So just remove the segment (by removing last point). + del pre_split_points[-1] + if len(pre_split_points) < 2: + # There are no segments before the intersection. + # So there's no pre-split polyline. + return None + + return pygplates.PolylineOnSphere(pre_split_points) # type: ignore + + def _post_split_polyline( + self, polyline, polyline_segment_index, polyline_intersection + ): + # Get the points *after* the intersection (including the intersection). + # This will have at least two points (needed for creating a polyline). + post_split_points = [polyline_intersection] + post_split_points.extend(polyline[polyline_segment_index + 1 :]) + + # Remove any zero length segments next to the intersected start of the polyline. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + while True: + post_split_first_segment = pygplates.GreatCircleArc( # type: ignore + post_split_points[0], + post_split_points[1], + ) + if not post_split_first_segment.is_zero_length(): + break + + # Remove the first segment (by removing first point). + del post_split_points[0] + if len(post_split_points) < 2: + # All segments after the intersection are zero length. + # So there's no post-split polyline. + return None + + # We've encountered a non-zero length segment. Let's shorten it such that it will no longer intersect the intersection point. + # Because we've already detected the intersection and don't want it detected in subsequent intersection tests. + if ( + self.INTERSECTION_THRESHOLD_RADIANS + < post_split_first_segment.get_arc_length() + ): + # First segment is *longer* than the intersection threshold, so shorten it by that amount. + # Replace first point with rotated point. + post_split_points[0] = ( + pygplates.FiniteRotation( # type: ignore + post_split_first_segment.get_rotation_axis(), + self.INTERSECTION_THRESHOLD_RADIANS, + ) + * post_split_points[0] + ) + else: + # First segment is *shorter* than the intersection threshold. But we know it's not a zero length segment, + # so it should be far enough away from the intersection point. So just remove the segment (by removing first point). + del post_split_points[0] + if len(post_split_points) < 2: + # There are no segments after the intersection. + # So there's no post-split polyline. + return None + + return pygplates.PolylineOnSphere(post_split_points) # type: ignore def reconstruct_points_by_topologies( - rotation_features_or_model, - topology_features, + plate_reconstruction, 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, + continent_mask_filepath_format=None, ): - """Reconstruct points using the topological polygons.""" + """Reconstruct points using topologies.""" - topology_reconstruction = _ReconstructByTopologiesImpl( - rotation_features_or_model, - topology_features, + topology_reconstruction = ReconstructByTopologies( + plate_reconstruction, reconstruction_begin_time, reconstruction_end_time, reconstruction_time_interval, points, point_begin_times=point_begin_times, point_end_times=point_end_times, - point_plate_ids=point_plate_ids, - detect_collisions=detect_collisions, + continent_mask_filepath_format=continent_mask_filepath_format, ) return topology_reconstruction.reconstruct() diff --git a/gplately/oceans.py b/gplately/oceans.py index 91cbb317..754f6fcc 100644 --- a/gplately/oceans.py +++ b/gplately/oceans.py @@ -30,10 +30,7 @@ from . import grids, tools from .lib.reconstruct_by_topologies import ( - _ContinentCollision, - _DefaultCollision, - _ReconstructByTopologiesImpl, - _ReconstructByTopologicalModelImpl, + ReconstructByTopologies, ) from .lib.reconstruct_continents import ReconstructContinents from .parallel import get_num_cpus @@ -107,13 +104,6 @@ class SeafloorGrid(object): Reconstruction by topologies involves determining which points are active and inactive (collided with a continent or subducted at a trench) for each reconstruction time step. This is done using a hidden object in :class:`PlateReconstruction` called ``ReconstructByTopologies``. - If an ocean point with a certain velocity on one plate ID transitions into another rigid plate ID at another timestep (with another velocity), - the velocity difference between both plates is calculated. The point may have subducted/collided with a continent if this velocity difference - is higher than a specified velocity threshold (which can be controlled with ``subduction_collision_parameters``). To ascertain whether the point - should be deactivated, a displacement test is conducted. If the proximity of the point's previous time position to the polygon boundary it is - approaching is higher than a set distance 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. Otherwise, it is deemed inactive and deleted from the ocean basin mesh. - With each reconstruction time step, points from mid-ocean ridges (which have more accurate spreading rates and attributed valid times) will spread across the ocean floor. Eventually, points will be pushed into continental boundaries or subduction zones, where they are deleted. Ideally, all initial ocean points (from the Stripy icosahedral mesh) should be deleted over time. However, not all will be deleted - such points typically reside near continental boundaries. @@ -196,6 +186,15 @@ def wrapper(self, *args, **kwargs): ) kwargs["continent_polygon_features"] = continent_polygon_features + # + # Handle the deprecated 'subduction_collision_parameters' argument. + # + if kwargs.pop("subduction_collision_parameters", None): + warnings.warn( + "`subduction_collision_parameters` keyword argument has been deprecated, it is no longer used", + DeprecationWarning, + ) + # # Handle the deprecated 'zval_names' argument. # @@ -223,7 +222,6 @@ def __init__( ridge_sampling: float = 0.5, extent: Tuple = (-180, 180, -90, 90), grid_spacing: float = 0.1, - subduction_collision_parameters=(5.0, 10.0), initial_ocean_mean_spreading_rate: float = 75.0, resume_from_checkpoints=False, continent_polygon_features=None, @@ -261,8 +259,6 @@ def __init__( grid_spacing : float, default=0.1 The degree spacing/interval with which to space grid points across all masking and final grids. If ``grid_spacing`` is provided, all grids will use it. If not, ``grid_spacing`` defaults to 0.1. - subduction_collision_parameters : len-2 tuple of float, default=(5.0, 10.0) - A 2-tuple of (threshold velocity delta in kms/my, threshold distance to boundary per My in kms/my) initial_ocean_mean_spreading_rate : float, default=75. A spreading rate to uniformly allocate to points that define the initial ocean basin. These points will have inaccurate ages, but most of them will be phased @@ -334,7 +330,6 @@ def __init__( # Topological parameters self.refinement_levels = refinement_levels self.ridge_sampling = ridge_sampling - self.subduction_collision_parameters = subduction_collision_parameters self.initial_ocean_mean_spreading_rate = initial_ocean_mean_spreading_rate # Gridding parameters @@ -1145,22 +1140,23 @@ def _calc_shifted_mor_point_and_spreading_rate( return shifted_mor_points, np.array(point_spreading_rates) + def reconstruct_by_topological_model(self): + """Alias for :meth:`reconstruct_by_topologies`. + + Introduced in version ``2.0``, this method used to use `pygplates.TopologicalModel`_ class to reconstruct seed points. + Now it's just an alias for :meth:`reconstruct_by_topologies` and is deprecated. + + .. deprecated:: 2.1 + + Use :meth:`reconstruct_by_topologies` instead. + """ + self.reconstruct_by_topologies() + def reconstruct_by_topologies(self): """Obtain all active ocean seed points which are points that have not been consumed at subduction zones or have not collided with continental polygons. Active points' latitudes, longitues, seafloor ages, spreading rates and all other general z-values are saved to a gridding input file (.npz). """ - self._reconstruct_by_topologies_impl(use_topological_model=False) - - def reconstruct_by_topological_model(self): - """Use `pygplates.TopologicalModel`_ class to reconstruct seed points. - This method is an alternative to :meth:`reconstruct_by_topologies()` which uses Python code to do the reconstruction. - - .. _pygplates.TopologicalModel: https://www.gplates.org/docs/pygplates/generated/pygplates.topologicalmodel - """ - self._reconstruct_by_topologies_impl(use_topological_model=True) - - def _reconstruct_by_topologies_impl(self, use_topological_model): logger.info("Preparing to reconstruct ocean seed points using topologies...") @@ -1211,15 +1207,12 @@ def _reconstruct_by_topologies_impl(self, use_topological_model): partial( _build_and_reconstruct_ocean_seed_points_parallel, seafloor_grid=self, - use_topological_model=use_topological_model, ), self._times, ) else: for time in self._times: - self._build_and_reconstruct_ocean_seed_points( - time, use_topological_model=use_topological_model - ) + self._build_and_reconstruct_ocean_seed_points(time) logger.info( "Aggregating reconstructed ocean seed point data for gridding input..." @@ -1305,9 +1298,7 @@ def _reconstruct_by_topologies_impl(self, use_topological_model): all_reconstructed_seed_point_data ) - def _build_and_reconstruct_ocean_seed_points( - self, from_time, use_topological_model - ): + def _build_and_reconstruct_ocean_seed_points(self, from_time): """Creates ocean seed points at `from_time` and reconstructs them to `min_time`. Ocean seed points can be either the *initial* ocean seed points at `from_time=max_time` or @@ -1377,48 +1368,18 @@ def _build_and_reconstruct_ocean_seed_points( # Begin the reconstruction by topology process. # - # Specify the default collision detection region as subduction zones - default_collision = _DefaultCollision( - feature_specific_collision_parameters=[ - ( - pygplates.FeatureType.gpml_subduction_zone, - self.subduction_collision_parameters, - ) - ] - ) - # In addition to the default subduction detection, also detect continental collisions - collision_spec = _ContinentCollision( - # This filename string should not have a time formatted into it - this is - # taken care of later. - self.continent_mask_filepath, - default_collision, - verbose=False, + # Call the reconstruct-by-topologies object. + topology_reconstruction = ReconstructByTopologies( + self.plate_reconstruction, + from_time, + self._min_time, + self._ridge_time_step, + points, + point_begin_times=appearance_times, + # This filename string should not have a time formatted into it - this is taken care of later... + continent_mask_filepath_format=self.continent_mask_filepath, ) - # Call the reconstruct by topologies object - if use_topological_model: - topology_reconstruction = _ReconstructByTopologicalModelImpl( - self.plate_reconstruction.rotation_model, - self.plate_reconstruction.topology_features, - from_time, - self._min_time, - self._ridge_time_step, - points, - point_begin_times=appearance_times, - detect_collisions=collision_spec, - ) - else: - topology_reconstruction = _ReconstructByTopologiesImpl( - self.plate_reconstruction.rotation_model, - self.plate_reconstruction.topology_features, - from_time, - self._min_time, - self._ridge_time_step, - points, - point_begin_times=appearance_times, - detect_collisions=collision_spec, - ) - # Initialise the reconstruction. topology_reconstruction.begin_reconstruction() @@ -1432,48 +1393,41 @@ def _build_and_reconstruct_ocean_seed_points( # This is the same size as 'points', which is either the initial ocean points at `time=max_time` or the MOR seed points at `time`. # Here 'current_points', despite being the same size, differs in that any *inactive* points at the current time # are None and the *active* points are at their reconstructed position at the current time. - current_points = topology_reconstruction.get_all_current_points() - - # Get the indices of the currently *active* points into all the points. - # - # Note: Points that are None represent currently *inactive* points. - # These are points that have either: - # 1) not been activated yet because the current time is older than their appearance time, or - # 2) have been deactivated because the current time is younger than their dissappearance time, or - # 3) have been deactivated through collision detection (ie, collided with a continent or trench). - # However, note that (2) does not apply here because the points currently don't have a hardwired disappearance time - # (it's implicitly '-inf'). Instead we rely on collision detection to deactivate points. - # And note that (1) does not apply here either, because the initial ocean seed points all exist at 'time=max_time' and - # all MOR seed points created at 'time' also exist at 'time'. - # - active_point_indices = np.fromiter( - ( - point_index - for point_index, point in enumerate(current_points) - if point is not None - ), - dtype=np.int32, - ) + active_points = topology_reconstruction.get_current_active_points() # logger.debug( - # f"At {time:.2f} Ma, {len(active_point_indices)} of the {len(points)} points created at {from_time:.2f} Ma are still active." + # f"At {time:.2f} Ma, {len(active_points)} of the {len(points)} points created at {from_time:.2f} Ma are still active." # ) # Store the reconstructed seed point data if any seed points are currently active. - if len(active_point_indices) > 0: + if active_points: # Latitudes and longitudes of the active points. # Store as 32-bit floating-point to save memory/disk-space. - active_latitudes = np.empty(len(active_point_indices), dtype=np.float32) - active_longitudes = np.empty( - len(active_point_indices), dtype=np.float32 + active_latitudes = np.empty(len(active_points), dtype=np.float32) + active_longitudes = np.empty(len(active_points), dtype=np.float32) + for index, active_point in enumerate(active_points): + active_latitudes[index], active_longitudes[index] = ( + active_point.to_lat_lon() + ) + + # Get the indices of the currently *active* points into all the points (active and inactive). + # + # Note: Points that are None represent currently *inactive* points. + # These are points that have either: + # 1) not been activated yet because the current time is older than their appearance time, or + # 2) have been deactivated because the current time is younger than their dissappearance time, or + # 3) have been deactivated through collision detection (ie, collided with a continent or trench). + # However, note that (2) does not apply here because the points currently don't have a hardwired disappearance time + # (it's implicitly '-inf'). Instead we rely on collision detection to deactivate points. + # And note that (1) does not really apply here either, because the initial ocean seed points all exist at the + # start time 'from_time=max_time' and all MOR seed points created at 'from_time' also exist at 'from_time'. + # + active_point_indices = ( + topology_reconstruction.get_current_active_point_indices() ) - for index, point_index in enumerate(active_point_indices): - active_latitudes[index], active_longitudes[index] = current_points[ - point_index - ].to_lat_lon() - # Store the active reconstructed seed point locations their active point indices for the current time. + # Store the active reconstructed seed point locations and their active point indices for the current time. # # Use a time "index" (associated with the current time - in the range `min_time <= time <= from_time`). time_index = topology_reconstruction.get_current_time_index() @@ -2012,9 +1966,6 @@ def _generate_debug_files_containing_reconstructed_ocean_seed_point_data_paralle def _build_and_reconstruct_ocean_seed_points_parallel( time, seafloor_grid, - use_topological_model, ): # Dispatch from parallel helper function to SeafloorGrid class method. - seafloor_grid._build_and_reconstruct_ocean_seed_points( - time, use_topological_model=use_topological_model - ) + seafloor_grid._build_and_reconstruct_ocean_seed_points(time)