iceberg sea level correction

README.md

@@ -42,21 +42,21 @@ Update when v.1.0 uploaded to Zendoo.
  The use of the pDEMtools package can be cited as follows:
 
 
-> Chudley, T. R. and Howat, I. M. (2024) pDEMtools (vX.X). GitHub. https://github.com/trchudley/pDEMtools
+> Chudley, T. R. and Howat, I. M. (2024) pDEMtools (vX.X.X). GitHub. https://github.com/trchudley/pDEMtools
 
 or by using `bibtex`:
 
 ```
 @software{pDEMtools
-   author = {Chudley, Thomas R. and Howat, Ian M.}, title = {pDEMtools}, year = 2024, publisher = {GitHub}, version = {X.X}, url = {https://github.com/trchudley/pDEMtools} 
+   author = {Chudley, Thomas R. and Howat, Ian M.}, title = {pDEMtools}, year = 2024, publisher = {GitHub}, version = {X.X.X}, url = {https://github.com/trchudley/pDEMtools} 
 }
 ```
 
 When using ArcticDEM and REMA products, please [cite](#refererences) the datasets appropriately and [acknowledge](#acknowledgements) the PGC.
 
 Several algorithms implemented in the library were developed by others. These will be highlighted in the documentation, and the original authors should be properly cited when used. For example:
 
-> We masked sea ice and melange following the method of Shiggins _et al._ (2023) as implemented in pDEMtools (Chudley, 2023).
+> We masked sea ice and melange following the method of Shiggins _et al._ (2023) as implemented in pDEMtools (Chudley and Howat, 2024).
 
 <!-- # To do
 

docs/getting_started/cite.md

@@ -2,18 +2,18 @@
 
 A software paper is being prepared for the [Journal of Open Source Software](https://joss.theoj.org/). In the meantime, the use of the pDEMtools package can be cited as follows:
 
-> Chudley, T. R. and Howat, I. M. (2024) pDEMtools (vX.X). GitHub. https://github.com/trchudley/pDEMtools
+> Chudley, T. R. and Howat, I. M. (2024) pDEMtools (vX.X.X). GitHub. https://github.com/trchudley/pDEMtools
 
 or by using `bibtex`:
 
 ```
 @software{pDEMtools
-   author = {Chudley, Thomas R. and Howat, Ian M.}, title = {pDEMtools}, year = 2024, publisher = {GitHub}, version = {X.X}, url = {https://github.com/trchudley/pDEMtools} 
+   author = {Chudley, Thomas R. and Howat, Ian M.}, title = {pDEMtools}, year = 2024, publisher = {GitHub}, version = {X.X.X}, url = {https://github.com/trchudley/pDEMtools} 
 }
 ```
 
 When using ArcticDEM and REMA products, please [cite](../appendix/references.md) the datasets appropriately and [acknowledge](../appendix/acknowledgements.md) the PGC.
 
 Several algorithms implemented in the library were developed by others. These will be highlighted in the documentation, and the original authors should be properly cited when used. For example:
 
-> We masked sea ice and melange following the method of Shiggins _et al._ (2023) as implemented in pDEMtools (Chudley, 2023).
+> We masked sea ice and melange following the method of Shiggins _et al._ (2023) as implemented in pDEMtools (Chudley and Howat, 2024).

docs/getting_started/supplementary_datasets.md

@@ -6,6 +6,6 @@ To make the most of pDEMtools, two supplementary datasets must available locally
 
 The first is the ArcticDEM or REMA strip index made available by the PGC  (in GeoParquet format), used by the `search` function. Strip index GeoParquet files can be downloaded from the PGC ([Greenland](https://data.pgc.umn.edu/elev/dem/setsm/ArcticDEM/indexes/), [Antarctica](https://data.pgc.umn.edu/elev/dem/setsm/REMA/indexes/)). To enable rapid searching, *please download the GeoParquet file format*: these files end in `_gpqt.zip`. Unzip them before use.
 
-## BedMachine v5
+## BedMachine
 
 The second is the Greenland BedMachine (v5; Morlighem _et al._ 2022a) or Antarctica BedMachine (v3; Morlighem _et al._ 2022b), which is the default geoid/bedrock mask used by the geoid correction and coregistration functions (NB: for applications outside of the ice sheets, it is possible to use your own geoid/bedrock mask). BedMachine can be downloaded from the NSIDC ([Greenland](https://nsidc.org/data/idbmg4/versions/5), [Antarctica](https://nsidc.org/data/nsidc-0756/versions/3)).

paper/paper.md

@@ -46,7 +46,7 @@ The `pdemtools` `search()` tool and `load` module allow for convenient access to
 
 After download, there exist a number of (pre-)processing steps that are near universally common in topographic analyses. These include geoid-correction, co-registration of time-series data, and/or the construction of terrain parameters such as hillshade, slope, aspect, and curvature. `pdemtools` contains pre-built functions to perform these processing steps, as well as further functionality specific to ArcticDEM and REMA use cases. For instance, we include functions to quickly extract the EIGEN-6C4 geoid [@foerste_geoid_2014] and Greenland/Antartctic bedrock masks directly from local versions of the Greenland and Antarctic BedMachine datasets [@morlighem_icebridge_2020; @morlighem_measures_2020], reprojecting and resampling the data to match the target DEM. Options for ingesting user-provided data are also provided. Additionally, partial deriviatives of the surface used to calculate terrain parameters ($\frac{\partial z}{\partial x}$, $\frac{\partial z}{\partial y}$, $\frac{\partial^2 z}{\partial x^2}$, $\frac{\partial^2 z}{\partial y^2}$, $\frac{\partial^2 z}{\partial x \partial y}$) are calculated following @florinsky_computation_2009, as opposed to more common methods such as @zevenbergen_quantitative_1987. The newer approach computes partial derivatives of elevation based on fitting a third-order polynomial, by the least-squares approach, to a 5 $\times$ 5 window as opposed to the more common 3 $\times$ 3 window. This is more appropriate for high-resolution DEMs: curvature over a 10 m window for the 2 m resolution ArcticDEM/REMA strips will lead to a local denoising effect that limits the impact of noise common in high-resolution photogrammetric products. These methods are also adapted into our co-registration routine, which otherwise follows the commonly used approach of @nuth_coregistration_2011.
 
-We aim to grow `pdemtools` by implementing new methods developed by the ArcticDEM and REMA research community. For instance, we currently include sea-level-filtering and iceberg detection routines outlined by @shiggins_automated_2023, and invite community contributions or requests of other routines that will be of use to users of `pdemtools`. Ongoing research projects making use of `pdemtools` are applying it to the mapping of crevasses, ice cliff heights, and subglacial lakes, as well is the initiation of ice sheet models. It has also been used within training exercises at the 2024 Polar Geospatial Center Data Workshop, contributing to a growing international network of `pdemtools` users.
+We aim to grow `pdemtools` by implementing new methods developed by the ArcticDEM and REMA research community. For instance, we currently include sea-level-filtering and iceberg detection routines outlined by @shiggins_automated_2023, and invite community contributions or requests of other routines that will be of use to users of `pdemtools`. Ongoing research projects making use of `pdemtools` are applying ArcticDEM and REMA data to the mapping of crevasses, ice cliff heights, and subglacial lakes, as well as the initiation of ice sheet models. It has also been used within training exercises at the 2024 Polar Geospatial Center Data Workshop, contributing to a growing international network of `pdemtools` users.
 
 
 # Documentation
@@ -56,7 +56,7 @@ Package documentation is available at the `pdemtools` [readthedocs](https://pdem
 
 # Acknowledgements
 
-Chudley is supported by a Leverhulme Early Career Fellowship (ECF-2022-589). ArcticDEM data are provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, and 1542736. REMA data are provided by the Byrd Polar and Climate Research Center and the Polar Geospatial Center under NSF-OPP awards 1543501, 1810976, 1542736, 1559691, 1043681, 1541332, 0753663, 1548562, 1238993 and NASA award NNX10AN61G. Computer time provided through a Blue Waters Innovation Initiative. DEMs produced using data from Maxar. We are grateful to the PGC for making index strips in the parquet format available from their website; to Connor Shiggins for his input on the iceberg detection function; to Andrew Sole and Adrien Wehrle for identifying and helping to resolve bugs in early versions of `pdemtools`; and to the many colleagues who helped to test this package in the early phases of development.
+Chudley is supported by a Leverhulme Early Career Fellowship (ECF-2022-589). ArcticDEM data are provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, and 1542736. REMA data are provided by the Byrd Polar and Climate Research Center and the Polar Geospatial Center under NSF-OPP awards 1543501, 1810976, 1542736, 1559691, 1043681, 1541332, 0753663, 1548562, 1238993 and NASA award NNX10AN61G. Computer time provided through a Blue Waters Innovation Initiative. DEMs produced using data from Maxar. We are grateful to the PGC for making index strips in the parquet format available from their website; to Connor Shiggins for his input on the iceberg detection function; to Andrew Sole and Adrien Wehrle for identifying and helping to resolve bugs in early versions of `pdemtools`; and to the colleagues who helped to test this package in the early phases of development.
 
 
 # References

src/pdemtools/_accessor.py

@@ -410,6 +410,7 @@ def mask_ocean(
         candidate_height_thresh_m: Optional[float] = 10,
         candidate_area_thresh_km2: Optional[float] = 1,
         near_sealevel_thresh_m: Optional[float] = 5,
+        return_sealevel_as_zero: Optional[bool] = False,
         return_mask: Optional[bool] = False,
     ) -> DataArray:
         """Returns a DEM with mélange/ocean regions filtered out, adapting the method
@@ -422,11 +423,15 @@ def mask_ocean(
             as SL, in m, defaults to 10
         :type candidate_height_thresh_m: float
         :param candidate_area_thresh_km2: Minimum area beneath `candidate_height_thresh_m`
-            to be considered for sea level assessment, in km^2, defaults to 1
+            to be considered for sea level assessment, in km^2. Defaults to 1
         :type candidate_area_thresh_km2: float
         :param near_sealevel_thresh_m: Filter out regions below this value, in metres above
-            sea level, defaults to 5
+            sea level. Defaults to 5
         :type near_sealevel_thresh_m: float
+        :param return_sealevel_as_zero: If True, subtracts that estimated sea level from
+            the DEM, returning new height values with the estimated sea level set to 0
+            m. Necessary for calculating accurate iceberg heights. Defaults to False.
+        :type return_sealevel_as_zero: bool
         :param return_mask: Return the sea level mask rather than the masked DEM.
             Defaults to False.
         :type return_mask: bool
@@ -457,13 +462,17 @@ def mask_ocean(
             if mask == None:
                 warn("No sea level detected in image. Returning `None` as mask object.")
             return mask
-
         else:
             if mask is None:
                 warn("No sea level detected in image. Returning original DEM")
                 return self._obj
             else:
-                return self._obj.where(mask)
+                if return_sealevel_as_zero == True:
+                    return self._obj.where(mask) - est_sea_level
+                elif return_sealevel_as_zero == False:
+                    return self._obj.where(mask)
+                else:
+                    raise ValueError("return_sealevel_as_zero must be True or False")
 
     def get_sea_level(
         self,
@@ -520,6 +529,13 @@ def mask_icebergs(
         identifying the size of connected groups of pixels and masking above/below a
         threshold.
 
+        By default, this function masks icebergs and retains terrestrial ice/land.
+        However, by setting `retain_icebergs = True`, the function will instead mask
+        ice/land and retain icebergs, allowing iceberg area and volume to be assessed.
+        For accurate volume assessment, the DEM 0 m value must be set to the estimated
+        sea level height from the `get_sea_level()` function.  This can be automated by
+        setting `return_sealevel_as_zero = True` in the mask_ocean() function.
+
         :param area_thresh_m2: Size threshold between icebergs and terrestrial ice/land, in
             m2. Defaults to 1e6 m2 (1 km2)
         :type area_thresh_m2: float | int