35 add dem array to terraincorrectionparameters

.gitignore

@@ -160,6 +160,9 @@ cython_debug/
 #.idea/
 gsolve.code-workspace
 docs/generated/*
+docs/_build/*
+docs/source/api/api/*
+docs/zz/*
 .vscode/settings.json
 .vscode/launch.json
 /*.ipynb

docs/source/Tutorial_terrain_correction.md

@@ -55,16 +55,24 @@ inner_zone_params = TerrainCorrectionParameters(
     water_density=1030.0,  # in kg/m3
     distance_mask_type="radial",
     dem_source=grid_dir / dem_file1,
-    clamp_elevation=True,
+    compute_topography=True,
+    compute_bathymetry=False,
+    site_height_field="dem_elevation",
     name="8m_dem",
 )
 
+
 outer_zone_params = TerrainCorrectionParameters(
     min_dist=2160.0,  # in meters
     max_dist=21900.0,  # in meters
+    terrain_density=2670.0,  # in kg/m3
+    water_density=1030.0,  # in kg/m3
+    distance_mask_type="radial",
     dem_source=grid_dir / dem_file2,
+    compute_topography=True,
+    compute_bathymetry=False,
+    site_height_field="dem_elevation",
     name="200m_dem",
-    clamp_elevation=False,
 )
 ```
 

docs/source/api/core_utils.rst

@@ -18,6 +18,15 @@ Functions
 =========
 .. currentmodule:: gsolve.core.utils
 
+Type Checking
+-------------
+.. autosummary::
+   :toctree: api/
+   is_filepath_like
+   is_datetime_array
+   is_in_literal
+
+
 Data Conversion
 ---------------
 .. autosummary::
@@ -52,7 +61,7 @@ Survey Loop Handling
    loops_from_gaps
 
 Other Utilities
--------------------
+---------------
 .. autosummary::
    :toctree: api/
 

docs/source/terrain_corrections.md

@@ -16,6 +16,9 @@ Terrain correction calculations in gSolve works by estimating the gravity effect
 
 The user imports the appropriate topography (or bathymetry) grid and observation location data.  Terrain and water densities, and minimum and maximum calculation distances are also required to be supplied.
 
+```{tip} gSolve can handle most grid formats readable by GDAL.  Tiff is well tested and preferred.
+```
+
 There is a choice of ```radial``` distance or ```rectangular``` distance.
 
 The algorithm trims the DEM to a box around the stations plus the ```max_dist```.  This reduces the number of prisms needed to be created and reduces memory.
@@ -40,11 +43,13 @@ To avoid duplicate calculations the following tips are useful.
 
 In ```TerrainCorrectionParameters``` there is an option ```compute_bathymetry``` or ```compute_topography```. For grids that contain both topography and bathymetry this can be used to set which part is computed.  The parameter ```sea_level_elevation``` can be used to control the boundary between topography and bathymetry (default = 0).
 
-The ```Sites``` object has a method ```sample_elevation()``` which can be used to extract the value of the DEM at the station location.  This can be useful if the DEM values are preferred to the surveyed positions or if you want to convert station ellipsoidal heights to orthometric (assuming the DEM is in orthometric heights).  A new column (```output_col=```) is created in the sites object for the interpolated DEM height which can be used in the computations through setting the ```site_height_field``` in the ```corrector.compute()```.
+The ```Sites``` object has a method ```sample_elevation()``` which can be used to extract the value of the DEM at the station location.  This can be useful if the DEM values are preferred to the surveyed positions or if you want to convert station ellipsoidal heights to orthometric (assuming the DEM is in orthometric heights).  A new column (```output_col=```) is created in the sites object for the interpolated DEM height and is used by setting the ``site_height_field`` in the ```TerraincCorrectionParameters```.
 
-gSolve will read most grid formats readable by GDAL.
+```{tip} If you are manipulating grids outside of gsolve they must be returned to gsolve as xarray.DataArray objects not numpy objects.
+```
 
-Densities are supplied in kg/m3.
+```{tip} Densities are supplied in kg/m3.
+```
 
 ## Terrain correction calculations over land
 
@@ -89,5 +94,3 @@ For calculation over the oceans the grids used at GNS are:
 
 * 2000m grid for terrain correction calculations between 21900 and 167000 m
   (derived from the resampled NIWA 250m dem)
-
-Note the source of the dems should be acknowledged in reports and papers.

examples/scripts/terrain_correction_example.py

@@ -16,7 +16,7 @@
 
 # -*- coding: utf-8 -*-
 """
-Created on Fri Mar  3 17:54:24 2023
+Created on Fri Mar  3 17:54:24 2023.
 
 @author: craigm
 """
@@ -64,7 +64,9 @@
     water_density=1030.0,  # in kg/m3
     distance_mask_type="radial",
     dem_source=grid_dir / dem_file1,
+    compute_topography=True,
     compute_bathymetry=False,
+    site_height_field="dem_elevation",
     name="8m_dem",
 )
 
@@ -75,16 +77,16 @@
     water_density=1030.0,  # in kg/m3
     distance_mask_type="radial",
     dem_source=grid_dir / dem_file2,
+    compute_topography=True,
     compute_bathymetry=False,
+    site_height_field="dem_elevation",
     name="200m_dem",
 )
 # Set up the Terrain corrector
 corrector = TerrainCorrector(params=[inner_zone_params, outer_zone_params])
 
 # Calculate the corrections
-results = corrector.compute(
-    points=sites, show_progress=True, site_height_field="dem_elevation"
-)
+results = corrector.compute(points=sites, show_progress=True)
 
 # Save results to file
 results.to_excel(

pyproject.toml

@@ -27,7 +27,8 @@ classifiers = [
 dependencies = [
     "boule==0.6.*",
     "harmonica<0.8",
-    "pandas[excel,plot]<3",
+    "openpyxl>=3.1.5",
+    "pandas[plot]<3",
     "pydata-sphinx-theme>=0.16.1",
     "pygtide>=0.8.2",
     "pyhardisp>=0.2.1",
@@ -41,7 +42,6 @@ dev = ["pytest>=9.0.3", "pytest-cov>=7.1.0", "pytest-datadir>=1.8.0", "ruff"]
 
 [project.urls]
 Repository = "https://github.com/GNS-Science/gsolve"
-# Homepage = "https://github.com/GNS-Science/gsolve"
 Documentation = "https://gns-science.github.io/gsolve/"
 Issues = "https://github.com/GNS-Science/gsolve/issues"
 Changelog = "https://github.com/GNS-Science/gsolve/blob/main/CHANGELOG.md"
@@ -116,14 +116,15 @@ ignore = ["ANN101", "ANN102", "ANN003", "D105", "DOC502", "DOC501"]
 preview = true
 
 [tool.ruff.lint.per-file-ignores]
-"tests/**" = ["ANN001", "ANN201", "ANN202"]
+"tests/**" = ["ANN001", "ANN201", "ANN202", "D"]
+"examples/scripts/**" = ["D100"]
 
 [tool.ruff.lint.pydocstyle]
 convention = "numpy"
 
 [tool.burocrata]
 notice = '''
-# GSolve - gravity processing software.
+# gSolve - gravity processing software.
 # Copyright (c) 2026 Earth Sciences New Zealand.
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by

src/gsolve/core/_typing.py

@@ -30,7 +30,6 @@
 from numpy.typing import ArrayLike, NDArray
 from pandas import DataFrame, DatetimeIndex, Index, Series, Timestamp
 from pandas.api.typing import NaTType
-from pylab import ndarray
 
 # from pandas.api.typing.aliases import TimedeltaConvertibleTypes
 
@@ -68,12 +67,12 @@
     NDArray, NDArray, NDArray, NDArray, NDArray, float | np.float64 | None, NDArray
 ]
 
-FilePath: TypeAlias = str | PathLike
+type FilePath = str | PathLike
 
 # The following type aliases are copied/adapted from pandas to ensure
 # function parameters are compatible with pandas methods they are passed to
 
-Renamer: TypeAlias = Union[Mapping[Any, Hashable], Callable[[Any], Hashable]]
+Renamer: TypeAlias = Mapping[Any, Hashable] | Callable[[Any], Hashable]
 
 
 DateTimeConvertibleTypes: TypeAlias = Union[
@@ -92,7 +91,7 @@
     int | float | str | datetime.date | np.datetime64 | pd.Timestamp
 )
 
-DatetimeArray: TypeAlias = list | tuple | ndarray | Series | Index | DatetimeIndex
+type DatetimeArray = list | tuple | Series | Index | DatetimeIndex | np.ndarray
 DatetimeScalarOrArray: TypeAlias = DatetimeScalar | DatetimeArray
 
 TimedeltaScalar: TypeAlias = str | int | float | pd.Timedelta | datetime.timedelta
@@ -107,14 +106,15 @@
 ArrayOrCoords: TypeAlias = DatasetOrArray | Sequence[ArrayLike]
 Points2D: TypeAlias = tuple[FloatArray, FloatArray]
 Points3D: TypeAlias = tuple[FloatArray, FloatArray, FloatArray]
-TCorrDistanceMaskType: TypeAlias = Literal["radial", "rectangular"]
+Points3DTrue = tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]
+type TCorrDistanceMaskType = Literal["radial", "rectangular"]
 
 
 # protocols for select Gsolve classes
 @runtime_checkable
 class SitesLike(Protocol):
     data: pd.DataFrame
 
-    def get_points(
-        self, xcol: str, ycol: str, zcol: str
-    ) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: ...
+    def get_points(self, xcol: str, ycol: str, zcol: str = "") -> Points3DTrue: ...
+
+    def get_site_ids(self) -> NDArray[np.str_]: ...

src/gsolve/core/excel_io.py

@@ -226,6 +226,7 @@ def write_excel_worksheet(
 
     writer_kwargs: dict[str, Any] = {
         # "engine": "openpyxl",  # "xlsxwriter", "openpyxl", "xlwt"
+        "engine": "openpyxl",
         "if_sheet_exists": if_sheet_exists,
         "mode": "w",
     }