diff --git a/gplately/grids.py b/gplately/grids.py index bcf24ec3..85983ca3 100644 --- a/gplately/grids.py +++ b/gplately/grids.py @@ -161,6 +161,41 @@ def _guess_data_variable_name(cdf: netCDF4.Dataset, x_name: str, y_name: str) -> return None +def _is_a_common_name_for_longitude(name: str) -> bool: + """Return True if the `name` parameter is a possible common name for longitude.""" + return name in ["lon", "lons", "longitude", "x", "east", "easting", "eastings"] + + +def _is_a_common_name_for_latitude(name: str) -> bool: + """Return True if the `name` parameter is a possible common name for latitude.""" + return name in ["lat", "lats", "latitude", "y", "north", "northing", "northings"] + + +def _find_extent_from_data(data) -> Union[Tuple[float, float, float, float], None]: + """Try to find the extent from data. Return None if data doesn't contain coordinates. + As of 2025-12-10, only support xarray.DataArray.""" + extent = None + lons = None + lats = None + try: + for name in data.coords: + if not lats and _is_a_common_name_for_latitude(name): + lats = data.coords[name] + elif not lons and _is_a_common_name_for_longitude(name): + lons = data.coords[name] + if lons is not None and lats is not None: + extent = ( + float(lons.min()), + float(lons.max()), + float(lats.min()), + float(lats.max()), + ) + except Exception as ex: + logger.debug(ex) + return None + return extent + + def read_netcdf_grid( filename, return_grids: bool = False, @@ -294,7 +329,9 @@ def find_label(keys, labels): if (unique_grid == [0, 1]).all(): cdf_grid = cdf_grid.astype(bool) - if realign: + # we realign the grid to -180/180 when the longitudes are from 0 to 360 + # this is a temporary fix. we need a more sophisticated solution. + if np.max(cdf_lon) > 180: # realign longitudes to -180/180 dateline cdf_grid_z, cdf_lon, cdf_lat = _realign_grid(cdf_grid, cdf_lon, cdf_lat) else: @@ -1564,6 +1601,7 @@ def __init__( The raster data, either as a file path (:class:`str`) or array data or a ``Raster`` object. If a ``Raster`` object is specified then all other arguments are ignored except ``plate_reconstruction`` which, if it is not ``None``, will override the plate reconstruction of the ``Raster`` object. + The data parameter accepts `numpy.ndarray`, `xarray.DataArray` or or any object that can be converted to a `numpy.ndarray`. plate_reconstruction : PlateReconstruction A :class:`PlateReconstruction` object to provide the following essential components for reconstructing points. @@ -1588,7 +1626,7 @@ def __init__( 2-tuple e.g. resample=(resX, resY). time : float, default: 0.0 - The geological time the time-dependant raster data. + The geological time of the time-dependant raster data. origin : {'lower', 'upper'}, optional When ``data`` is an array, use this parameter to specify the origin @@ -1645,6 +1683,7 @@ def __init__( + "'{}'".format(key) ) + # if the "data" parameter is a "Raster" object if isinstance(data, self.__class__): self._data = data._data.copy() # Use specified plate reconstruction (if specified), @@ -1661,15 +1700,20 @@ def __init__( self.plate_reconstruction = plate_reconstruction - if time < 0.0: - raise ValueError("Invalid time: {}".format(time)) - time = float(time) - self._time = time + # get the geological time parameter for the time-dependant raster data + try: + time = float(time) + if time < 0.0: + raise ValueError() + self._time = time + except ValueError: + raise ValueError(f"Invalid time parameter: {time}") if data is None: raise TypeError("`data` argument (or `filename` or `array`) is required") + + # if the user has passed a NetCDF file path if isinstance(data, str): - # Filename self._filename = data self._data, lons, lats = read_netcdf_grid( data, @@ -1683,17 +1727,30 @@ def __init__( ) self._lons = lons self._lats = lats - else: - # numpy array + # if the "data" parameter is a numpy array or xarray.DataArray object self._filename = None extent = _parse_extent_origin(extent, origin) + + # try to extract the extent from input data + # if the extent from data is different from the extent parameter, use the extent from data + extent_from_data = _find_extent_from_data(data) + if extent_from_data is not None and extent != extent_from_data: + extent = extent_from_data + logger.info( + f"Raster.__init__(): Use the extent extracted from data: {extent}." + ) + data = _check_grid(data) self._data = np.array(data) self._lons = np.linspace(extent[0], extent[1], self.data.shape[1]) self._lats = np.linspace(extent[2], extent[3], self.data.shape[0]) - if realign: - # realign to -180,180 and flip grid + + # we realign the grid to -180/180 when the longitudes are from 0 to 360 + # this is a temporary fix. we need a more sophisticated solution. + # for example, some people may use (-360-0) or some other ranges for longitudes. It is unlikely, but possible. + if np.max(self._lons) > 180: + # realign to -180,180 and flip grid if needed self._data, self._lons, self._lats = _realign_grid( self._data, self._lons, self._lats ) diff --git a/tests-dir/debug_utils/369-raster-reconstruct-in-0-360-format.ipynb b/tests-dir/debug_utils/369-raster-reconstruct-in-0-360-format.ipynb new file mode 100644 index 00000000..a554e879 --- /dev/null +++ b/tests-dir/debug_utils/369-raster-reconstruct-in-0-360-format.ipynb @@ -0,0 +1,322 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "23da50d7-50fc-4b01-92f7-9c398afac849", + "metadata": {}, + "source": [ + "# Raster reconstruct with 0-360 data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "815cca41-de45-4b53-a313-9bff81b1255a", + "metadata": {}, + "outputs": [], + "source": [ + "import gplately\n", + "import pygmt\n", + "import numpy as np\n", + "import xarray as xr\n", + "import os\n", + "import plate_model_manager\n", + "\n", + "def gplately_raster_to_xarray(raster, raster_name, raster_units):\n", + "\n", + " # check the order of lats\n", + " if raster.lats.min() < raster.lats.max():\n", + " ds_raster = xr.Dataset(\n", + " data_vars=dict(\n", + " z=(\n", + " [\"lat\", \"lon\"],\n", + " raster.data,\n", + " {\n", + " \"long_name\": raster_name,\n", + " \"units\": raster_units,\n", + " \"actual_range\": np.array(\n", + " [np.nanmin(raster.data), np.nanmax(raster.data)],\n", + " dtype=np.float32,\n", + " ),\n", + " },\n", + " )\n", + " ),\n", + " coords=dict(\n", + " lon=(\n", + " [\"lon\"],\n", + " raster.lons.astype(\"float32\"),\n", + " {\n", + " \"long_name\": \"longitude\",\n", + " \"units\": \"degrees_east\",\n", + " \"standard_name\": \"longitude\",\n", + " \"actual_range\": np.array(\n", + " [raster.lons.min(), raster.lons.max()], dtype=np.float32\n", + " ),\n", + " },\n", + " ),\n", + " lat=(\n", + " [\"lat\"],\n", + " (raster.lats).astype(\"float32\"),\n", + " {\n", + " \"long_name\": \"latitude\",\n", + " \"units\": \"degrees_north\",\n", + " \"standard_name\": \"latitude\",\n", + " \"actual_range\": np.array(\n", + " [raster.lats.min(), raster.lats.max()], dtype=np.float32\n", + " ),\n", + " },\n", + " ),\n", + " ),\n", + " )\n", + " else:\n", + " # otherwise the lats are in DESCENDING order, which means we need to flip the data and lats (otherwise the grid will be upside down later)\n", + " ds_raster = xr.Dataset(\n", + " data_vars=dict(\n", + " z=(\n", + " [\"lat\", \"lon\"],\n", + " np.flipud(raster.data),\n", + " {\n", + " \"long_name\": raster_name,\n", + " \"units\": raster_units,\n", + " \"actual_range\": np.array(\n", + " [np.nanmin(raster.data), np.nanmax(raster.data)],\n", + " dtype=np.float32,\n", + " ),\n", + " },\n", + " )\n", + " ),\n", + " coords=dict(\n", + " lon=(\n", + " [\"lon\"],\n", + " raster.lons.astype(\"float32\"),\n", + " {\n", + " \"long_name\": \"longitude\",\n", + " \"units\": \"degrees_east\",\n", + " \"standard_name\": \"longitude\",\n", + " \"actual_range\": np.array(\n", + " [raster.lons.min(), raster.lons.max()], dtype=np.float32\n", + " ),\n", + " },\n", + " ),\n", + " lat=(\n", + " [\"lat\"],\n", + " np.flipud(raster.lats).astype(\"float32\"),\n", + " {\n", + " \"long_name\": \"latitude\",\n", + " \"units\": \"degrees_north\",\n", + " \"standard_name\": \"latitude\",\n", + " \"actual_range\": np.array(\n", + " [raster.lats.min(), raster.lats.max()], dtype=np.float32\n", + " ),\n", + " },\n", + " ),\n", + " ),\n", + " )\n", + "\n", + " return ds_raster" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "654df182-4673-4812-9030-a6730e974ddb", + "metadata": {}, + "outputs": [], + "source": [ + "plate_models = \"plate_models\"\n", + "if not os.path.exists(plate_models):\n", + " os.makedirs(plate_models)\n", + "\n", + "# import target plate model: Zahirovic et al. (2022).\n", + "try:\n", + " S12_model_pmm = plate_model_manager.PlateModelManager().get_model(\n", + " \"Seton2012\", data_dir=plate_models\n", + " )\n", + "except:\n", + " # if unable to connect to the servers, try to use the local files\n", + " S12_model_pmm = plate_model_manager.PlateModel(\n", + " model_name=\"Seton2012\", data_dir=plate_models, readonly=True\n", + " )\n", + "\n", + "# get the plate model components\n", + "S12_rotation_filename = S12_model_pmm.get_rotation_model()\n", + "S12_static_polygons_filename = S12_model_pmm.get_static_polygons()\n", + "S12_coastline_filename = S12_model_pmm.get_coastlines()\n", + "S12_topology_features = S12_model_pmm.get_topologies()\n", + "\n", + "# create a GPlately 'plate reconstruction' object\n", + "S12_model = gplately.PlateReconstruction(\n", + " S12_rotation_filename,\n", + " S12_topology_features,\n", + " S12_static_polygons_filename,\n", + " anchor_plate_id=0,\n", + ")\n", + "\n", + "seton_2012_age_grid = gplately.Raster(\n", + " data=S12_model_pmm.get_raster(\"AgeGrids\", 50), plate_reconstruction=S12_model\n", + ")\n", + "\n", + "# check longitude extent\n", + "seton_2012_age_grid.lons\n", + "\n", + "# convert to dataset\n", + "ds_seton_2012_age_grid = gplately_raster_to_xarray(seton_2012_age_grid, \"age\", \"myr\")\n", + "\n", + "#ds_seton_2012_age_grid.z.plot()\n", + "\n", + "# reformat longitudes into 0-360 format\n", + "ds_seton_2012_age_grid[\"lon\"] = (ds_seton_2012_age_grid.coords[\"lon\"]) % 360\n", + "ds_seton_2012_age_grid = ds_seton_2012_age_grid.sortby(ds_seton_2012_age_grid.lon)\n", + "ds_seton_2012_age_grid.z.plot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ca1fad3-a654-4b08-84ed-343e9151a62a", + "metadata": {}, + "outputs": [], + "source": [ + "# create a gplately plotting object with our plate model\n", + "gplot = gplately.PlotTopologies(\n", + " S12_model,\n", + " coastlines=S12_coastline_filename,\n", + " continents=S12_static_polygons_filename,\n", + " COBs=None,\n", + " time=0,\n", + ")\n", + "gplot.time = 50\n", + "\n", + "# ---- actually plot\n", + "region = \"d\"\n", + "width = 8\n", + "projection = \"N180/%sc\" % width\n", + "\n", + "pen = \"1p\"\n", + "alpha = 25\n", + "\n", + "label_font = \"10p,Helvetica,black\"\n", + "label_offset = \"j0/0c\"\n", + "label_position = \"TL\"\n", + "\n", + "fig = pygmt.Figure()\n", + "\n", + "pygmt.makecpt(cmap=\"magma\", reverse=True, series=[0, 150, 1])\n", + "\n", + "pygmt.config(\n", + " FONT_ANNOT=6.5,\n", + " FONT_LABEL=8,\n", + " MAP_TICK_PEN=\"0.75p\",\n", + " MAP_FRAME_PEN=\"0.75p\",\n", + " MAP_TICK_LENGTH_PRIMARY=\"3p\",\n", + ")\n", + "fig.basemap(projection=projection, region=region, frame=\"lrtb\")\n", + "fig.grdimage(grid=ds_seton_2012_age_grid.z)\n", + "fig.plot(\n", + " gplot.get_continents(), pen=\"0.05p,black\", region=region, projection=projection\n", + ")\n", + "fig.text(\n", + " text=\"%s Ma\" % 50,\n", + " position=label_position,\n", + " no_clip=True,\n", + " font=label_font,\n", + " offset=label_offset,\n", + ")\n", + "fig.colorbar(frame=[\"x+lAge (myr)\"], position=\"JRM+o0.5c/0c+w%sc/0.2c+ef\" % (width - 4))\n", + "fig.show(width=1000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f3fa924d-e30f-4693-a9d4-e05892db5237", + "metadata": {}, + "outputs": [], + "source": [ + "# Reconstruct to present day\n", + "\n", + "# after PR #372 https://github.com/GPlates/gplately/pull/372\n", + "# no need to specify \"extent=(0, 360, -90, 90)\" \n", + "# the extent will be extracted from the xarray.DataArray object\n", + "raster_50 = gplately.Raster(data=ds_seton_2012_age_grid.z, time=50)\n", + "raster_50.plate_reconstruction = S12_model\n", + "\n", + "raster_0 = raster_50.reconstruct(0) # reconstruct to present-day\n", + "\n", + "ds_raster_0 = gplately_raster_to_xarray(raster_0, \"age\", \"myr\")\n", + "\n", + "gplot = gplately.PlotTopologies(\n", + " S12_model,\n", + " coastlines=S12_coastline_filename,\n", + " continents=S12_static_polygons_filename,\n", + " COBs=None,\n", + " time=0,\n", + ")\n", + "gplot.time = 0\n", + "\n", + "# ---- actually plot\n", + "region = \"d\"\n", + "width = 8\n", + "projection = \"N180/%sc\" % width\n", + "label_font = \"10p,Helvetica,black\"\n", + "label_offset = \"j0/0c\"\n", + "label_position = \"TL\"\n", + "\n", + "fig = pygmt.Figure()\n", + "\n", + "pygmt.makecpt(cmap=\"magma\", reverse=True, series=[0, 150, 1])\n", + "\n", + "pygmt.config(\n", + " FONT_ANNOT=6.5,\n", + " FONT_LABEL=8,\n", + " MAP_TICK_PEN=\"0.75p\",\n", + " MAP_FRAME_PEN=\"0.75p\",\n", + " MAP_TICK_LENGTH_PRIMARY=\"3p\",\n", + ")\n", + "fig.basemap(projection=projection, region=region, frame=\"lrtb\")\n", + "fig.grdimage(grid=ds_raster_0.z)\n", + "fig.plot(\n", + " gplot.get_continents(), pen=\"0.05p,black\", region=region, projection=projection\n", + ")\n", + "fig.text(\n", + " text=\"%s Ma\" % 0,\n", + " position=label_position,\n", + " no_clip=True,\n", + " font=label_font,\n", + " offset=label_offset,\n", + ")\n", + "fig.colorbar(frame=[\"x+lAge (myr)\"], position=\"JRM+o0.5c/0c+w%sc/0.2c+ef\" % (width - 4))\n", + "fig.show(width=1000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffbbceab-f0ac-46f5-af0e-e91073df6143", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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/tests-dir/debug_utils/github-issue-287-and-289.ipynb b/tests-dir/debug_utils/github-issue-287-and-289.ipynb index 6130d6b9..b7f6f6db 100644 --- a/tests-dir/debug_utils/github-issue-287-and-289.ipynb +++ b/tests-dir/debug_utils/github-issue-287-and-289.ipynb @@ -18,6 +18,9 @@ "from plate_model_manager import PlateModelManager, PlateModel\n", "import xarray as xr\n", "\n", + "# https://github.com/GPlates/gplately/issues/287\n", + "# https://github.com/GPlates/gplately/issues/289\n", + "\n", "workdir = \"debug-folder-github-issue-287\"\n", "Path(workdir).mkdir(parents=True, exist_ok=True)\n", "\n", @@ -49,9 +52,11 @@ " static_polygons=plate_model.get_layer(\"StaticPolygons\"),\n", ")\n", "\n", + "# After PR #372 https://github.com/GPlates/gplately/pull/372\n", + "# the \"realign=True\" parameter is not needed anymore\n", "raster = gplately.Raster(\n", - " data=miocene_file, time=15, plate_reconstruction=r_model, realign=True\n", - ")\n", + " data=miocene_file, time=15, plate_reconstruction=r_model)\n", + "\n", "reconstructed_raster: gplately.Raster = raster.reconstruct(\n", " 0, threads=4, partitioning_features=plate_model.get_layer(\"StaticPolygons\")\n", ") # type: ignore\n", @@ -95,7 +100,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.9" + "version": "3.11.11" } }, "nbformat": 4,