Feature:gia2plmt-SELEN compatibility

aux/addanchors.m

@@ -0,0 +1,116 @@
+%% ADDANCHORS
+% Add anchors to horizontal or vertical segments of a polygon to ensure that the segments are not too long.
+% This is useful when trying to plot the polygon on a map, as the segments are not displayed as great circle arcs.
+%
+% Syntax
+%   p = addanchors(p)
+%   lonlat = addanchors(lonlat)
+%   [lon, lat] = addanchors(lon, lat)
+%
+% Input arguments
+%   p - A polyshape object
+%   lonlat - A 2-column matrix of longitudes and latitudes
+%   lon, lat - Vectors of longitudes and latitudes
+%
+% Output arguments
+%   p - A polyshape object
+%   lonlat - A 2-column matrix of longitudes and latitudes
+%   lon, lat - Vectors of longitudes and latitudes
+%
+% Notes
+%   This function was originally part of the slepian_ulmo package.
+%
+% Last modified by
+%   2024/10/15, En-Chi Lee (williameclee@arizona.edu)
+
+function varargout = addanchors(varargin)
+    isPoly = false;
+
+    if nargin == 1
+
+        if isa(varargin{1}, 'polyshape')
+            isPoly = true;
+            varargin{1} = varargin{1}.Vertices;
+            varargin{1} = closecoastline(varargin{1});
+        elseif size(varargin{1}, 1) == 2
+            varargin{1} = varargin{1}';
+        elseif size(varargin{1}, 2) ~= 2
+            error('Invalid input argument')
+        end
+
+        lon = varargin{1}(:, 1);
+        lat = varargin{1}(:, 2);
+    elseif nargin == 2
+        lon = varargin{1};
+        lat = varargin{2};
+    else
+        error('Invalid number of input arguments')
+    end
+
+    lon = lon(:);
+    lat = lat(:);
+
+    maxSep = 5;
+    tol = 1;
+
+    diffLat = abs(diff(lat));
+    diffLon = abs(diff(lon));
+
+    needsAnchor = find( ...
+        (diffLat > maxSep & diffLon <= tol) | ...
+        (diffLon > maxSep & diffLat <= tol));
+    difLength = max(diffLat(needsAnchor), diffLon(needsAnchor));
+
+    lata = lat;
+    lona = lon;
+
+    if isempty(needsAnchor)
+
+        if nargout == 1
+
+            if isPoly
+                varargout = {p};
+            else
+                varargout = {[lon, lat]};
+            end
+
+        else
+            varargout = {lon, lat};
+        end
+
+        return
+    end
+
+    for i = length(needsAnchor):-1:1
+        latInterp = interp1([0, 1], ...
+            [lata(needsAnchor(i)), lata(needsAnchor(i) + 1)], ...
+            linspace(0, 1, round(difLength(i) / tol)));
+        lata = [lata(1:needsAnchor(i) - 1); latInterp'; lata(needsAnchor(i) + 1:end)];
+        lonInterp = interp1([0, 1], ...
+            [lona(needsAnchor(i)), lona(needsAnchor(i) + 1)], ...
+            linspace(0, 1, round(difLength(i) / tol)));
+        lona = [lona(1:needsAnchor(i) - 1); lonInterp'; lona(needsAnchor(i) + 1:end)];
+    end
+
+    if nargout == 0
+        figName = sprintf('Polygon with more anchors (%s)', upper(mfilename));
+        figure(998)
+        set(gcf, 'Name', figName, 'NumberTitle', 'off')
+        clf
+        plotqdm(lona, lata)
+        hold on
+        plot(lona, lata, 'k.-')
+        hold off
+    elseif nargout == 1
+
+        if isPoly
+            varargout = {polyshape(lona, lata, "Simplify", false)};
+        else
+            varargout = {[lona, lata]};
+        end
+
+    else
+        varargout = {lona, lata};
+    end
+
+end

domains/icesheetPoly.m

@@ -0,0 +1,55 @@
+%% ICESHEETPOLY - Gets an approximate ice sheet polygon  
+%
+% Syntax
+%   [icePoly, iceLonlat] = ICESHEETPOLY(icesheet)
+%
+% Input arguments
+%   icesheet - Name of the ice sheet
+%       The following ice sheets are recognised (case-insensitive):
+%       - "Fennoscandian", "Fennoscandian Ice Sheet", "FIS"
+%       - "Laurentide", "Laurentide Ice Sheet", "LIS"
+%       - "Greenland", "Greenland Ice Sheet", "GIS"
+%       - "Patagonian", "Patagonian Ice Sheet", "PIS"
+%       - "Antarctic", "Antarctic Ice Sheet", "AIS"
+%       - "West Antarctic", "West Antarctic Ice Sheet", "WAIS"
+%       - "East Antarctic", "East Antarctic Ice Sheet", "EAIS"
+%
+% Output arguments
+%   icePoly - Polyshape object of the ice sheet
+%   iceLonlat - Longitude-latitude vertices of the ice sheet polygon
+%
+% Author
+%   2026/04/07, En-Chi Lee (williameclee@arizona.edu)
+
+function [icePoly, iceLonlat] = icesheetPoly(icesheet)
+
+    arguments (Input)
+        icesheet {mustBeTextScalar}
+    end
+
+    arguments (Output)
+        icePoly (1, 1) polyshape
+        iceLonlat (:, 2) double
+    end
+
+    switch lower(icesheet)
+        case {"fennoscandian", "fennoscandian ice sheet", "fis"}
+            iceLonlat = [0, 45; 30, 45; 95, 70; 95, 90; 0, 90; NaN, NaN; 360, 68; 344, 58; 344, 47; 360, 47];
+        case {"laurentide", "laurentide ice sheet", "lis"}
+            iceLonlat = [240, 90; 185, 70; 185, 50; 215, 50; 230, 30; 290, 30; 315, 45; 315, 55; 300, 60; 300, 70; 286, 77; 300, 90];
+        case {"greenland", "greenland ice sheet", "gis"}
+            iceLonlat = [300, 60; 300, 70; 286, 77; 300, 90; 355, 90; 355, 75; 340, 70; 325, 60; 315, 55];
+        case {"patagonian", "patagonian ice sheet", "pis"}
+            iceLonlat = [275, -35; 275, -60; 300, -60; 300, -35];
+        case {"antarctic", "antarctic ice sheet", "ais"}
+            iceLonlat = [0, -90; 360, -90; 360, -60; 0, -60];
+        case {"west antarctic", "west antarctic ice sheet", "wais"}
+            iceLonlat = [185, -60; 166, -78; 153, -83; 162, -85; 200, -87; 240, -85; 290, -85; 306, -83; 318, -83; 340, -74; 340, -60];
+        case {"east antarctic", "east antarctic ice sheet", "eais"}
+            iceLonlat = [0, -90; 0, -60; 185, -60; 166, -78; 153, -83; 162, -85; 200, -87; 240, -85; 290, -85; 306, -83; 318, -83; 340, -74; 340, -60; 360, -60; 360, -90];
+        otherwise
+            error("Unrecognised ice sheet name: %s", icesheet)
+    end
+
+    icePoly = polyshape(iceLonlat(:, 1), iceLonlat(:, 2));
+end

externalio/selen/convertSelenModel.m

@@ -0,0 +1,72 @@
+%% CONVERTSELENMODEL - Converts Selen outputs to the format used by SLEPIAN
+% Converts the geoid and vertical displacement outputs from the SELEN model
+% to the format used by SLEPIAN_DELTA. This is not intended to be called
+% directly by users.
+%
+% Syntax
+%   convertSelenModel(modelFolder)
+%
+% Input arguments
+%   modelFolder - folder containing the SELEN model outputs
+%
+% References
+%   Spada, Giorgio & Melini, Daniele. (2019). SELEN4 (SELEN version 4.0): a
+%       Fortran program for solving the gravitationally and topographically
+%       self-consistent sea-level equation in glacial isostatic adjustment
+%       modeling. Geosci. Model Dev. 
+%       https://doi.org/10.5194/gmd-12-5055-2019
+%       https://github.com/geodynamics/selen
+%   Tamisiea, Mark E. (2011). Ongoing glacial isostatic contributions to
+%       observations of sea level change. Geophys. J. Int.
+%       https://doi.org/10.1111/j.1365-246X.2011.05116.x
+%
+% Author
+%	2026/04/01, En-Chi Lee (williameclee@arizona.edu)
+
+function convertSelenModel(modelFolder)
+
+    arguments (Input)
+        modelFolder {mustBeFolder}
+    end
+
+    dataFolder = fullfile(modelFolder, "FPR");
+    % Detect the degree of the model from the model folder name
+    [~, modelFolderSelf] = fileparts(modelFolder);
+    modelName = regexp(modelFolderSelf, "RUN_([A-Za-z0-9-_]+)", "tokens", "once");
+    modelInfo = regexp(modelFolderSelf, "RUN_([A-Za-z0-9]+)-R(\d+)-L(\d+)-I(\d+)", "tokens", "once");
+
+    % Prepare grid for interpolation
+    L = str2double(modelInfo{3});
+    h = 360 / (2 * L);
+    lon = 0:h:360;
+    lat = -90:h:90;
+    [lonn, latt] = meshgrid(lon, lat);
+    %% Geoid and surface density
+    geoidPath = fullfile(dataFolder, "gdot.pix");
+    geoidXyz = readmatrix(geoidPath, "FileType", "text");
+    geoidXyz(:, 3) = geoidXyz(:, 3) / 1e3; % mm -> m
+    F = scatteredInterpolant(geoidXyz(:, 1), geoidXyz(:, 2), geoidXyz(:, 3), "nearest", "nearest");
+    geoidMesh = F(mod(lonn, 360), latt);
+    geoidPlm = xyz2plm_new(flip(geoidMesh), L);
+
+    geoidPlm(5, 3:4) = geoidPlm(5, 3:4) / 2.06; % See Appendix of Tamisiea (2011)
+
+    sdPlm = convertgravity(geoidPlm, "POT", "SD");
+
+    % Save the converted model
+    lmcosiM = sdPlm;
+    save(fullfile(modelFolder, sprintf("%s_SD.mat", modelName)), "lmcosiM");
+    lmcosiM = geoidPlm;
+    save(fullfile(modelFolder, sprintf("%s_POT.mat", modelName)), "lmcosiM");
+
+    %% Vertical displacement
+    vlmPath = fullfile(dataFolder, "udot.pix");
+    vlmXyz = readmatrix(vlmPath, "FileType", "text");
+    F = scatteredInterpolant(vlmXyz(:, 1), vlmXyz(:, 2), vlmXyz(:, 3), "nearest", "nearest");
+    vlmMesh = F(mod(lonn, 360), latt);
+    vlmPlm = xyz2plm_new(flip(vlmMesh), L);
+
+    % Save the converted model
+    lmcosiM = vlmPlm;
+    save(fullfile(modelFolder, sprintf("%s_VLM.mat", modelName)), "lmcosiM");
+end

externalio/selen/modSelenIceModel.m

@@ -0,0 +1,141 @@
+%% MODSELENICEMODEL - Modifies a ice model in SELEN format
+% This function modifies the ice model used in a SELEN run by changing the
+% volume change and time delay of a specified ice sheet. The modified model
+% is saved as a new file in the same format as the original model. This
+% function is intended to be used for sensitivity tests of the ice model in
+% SELEN.
+%
+% Syntax
+%   modSelenIceModel(model, icesheet, newModel, volumeChangeFrac, timeDelay)
+%
+% Input arguments
+%   model - Name of the SELEN model to modify (without the .pix extension)
+%       The model must be located in the DATA folder of the SELEN
+%       installation.
+%   icesheet (optional) - Name of the ice sheet to modify
+%       See ICESHEETPOLY for recognised ice sheet names.
+%       The default ice sheet to modify is "fennoscandian".
+%   newModel (optional) - Name of the new model to save
+%       If empty, the new name is programmatically generated based on the
+%       original model name and the modifications made.
+%   volumeChangeFrac (optional) - Fractional change in the volume of the
+%       specified ice sheet
+%       For example, a value of 0.1 means a 10% increase in volume, while a
+%       value of -0.2 means a 20% decrease in volume.
+%       The default value is 0, meaning no change in volume.
+%   timeDelay (optional) - Time delay in years for the specified ice sheet
+%       A positive value means the ice sheet changes later than in the
+%       original model, while a negative value means it changes earlier.
+%       The default value is 0, meaning no time delay.
+%
+% See also
+%   ICESHEETPOLY
+%
+% References
+%   Spada, Giorgio & Melini, Daniele. (2019). SELEN4 (SELEN version 4.0): a
+%       Fortran program for solving the gravitationally and topographically
+%       self-consistent sea-level equation in glacial isostatic adjustment
+%       modeling. Geosci. Model Dev.
+%       https://doi.org/10.5194/gmd-12-5055-2019
+%
+% Author
+%   2026/04/07, En-Chi Lee (williameclee@arizona.edu)
+
+function modSelenIceModel(model, icesheet, newModel, volumeChangeFrac, timeDelay)
+
+    arguments (Input)
+        model {mustBeTextScalar} = "i6g-R44r-i"
+        icesheet {mustBeTextScalar} = "fennoscandian"
+        newModel {mustBeTextScalar} = ""
+        volumeChangeFrac (1, 1) double {mustBeReal, mustBeFinite} = 0
+        timeDelay (1, 1) double {mustBeInteger, mustBeReal, mustBeFinite} = 0
+    end
+
+    modelPath = fullfile(getenv("SELEN"), "DATA", sprintf("%s.pix", model));
+    icesheet = icesheetCode(icesheet);
+
+    ice = readmatrix(modelPath, "FileType", "text");
+    isIcesheet = labelIceSheet(ice(:, 2:3), icesheet);
+
+    iceDiff = ice(:, 5:end) - ice(:, end);
+
+    if volumeChangeFrac ~= 0
+        iceDiff(isIcesheet, :) = iceDiff(isIcesheet, :) * (1 + volumeChangeFrac);
+    end
+
+    if timeDelay > 0
+        iceDiff(isIcesheet, timeDelay + 1:end) = iceDiff(isIcesheet, 1:end - timeDelay);
+
+        iceDiff(isIcesheet, 1:timeDelay) = repmat(iceDiff(isIcesheet, 1), [1, timeDelay]);
+
+    elseif timeDelay < 0
+        iceDiff(isIcesheet, 1:end + timeDelay) = iceDiff(isIcesheet, -timeDelay + 1:end);
+    end
+
+    iceNew = [ice(:, 1:4), ice(:, end) + iceDiff];
+
+    if ~isempty(newModel)
+        newModel = sprintf("%s-%s_V%s_T%s", newModel, upper(icesheet), ...
+            replace(replace(sprintf("%+03.0f", volumeChangeFrac * 100), "+", "p"), "-", "n"), ...
+            replace(replace(sprintf("%+03d", timeDelay), "+", "p"), "-", "n"));
+    end
+
+    newModelPath = fullfile(getenv("SELEN"), "DATA", newModel);
+
+    fid = fopen(newModelPath, "w");
+    cleanUp = onCleanup(@() fclose(fid));
+
+    for i = 1:size(iceNew, 1)
+        fprintf(fid, "%12d %20f %20f %25f     ", iceNew(i, 1:4));
+        fprintf(fid, " %11.0f", iceNew(i, 5:end));
+        fprintf(fid, "\n");
+    end
+
+end
+
+%% Subfunctions
+function code = icesheetCode(name)
+
+    arguments (Input)
+        name {mustBeTextScalar}
+    end
+
+    arguments (Output)
+        code (1, 1) string
+    end
+
+    switch lower(name)
+        case {"fennoscandian", "fennoscandian ice sheet", "fis"}
+            code = "FIS";
+        case {"laurentide", "laurentide ice sheet", "lis"}
+            code = "LIS";
+        case {"greenland", "greenland ice sheet", "gis"}
+            code = "GIS";
+        case {"patagonian", "patagonian ice sheet", "pis"}
+            code = "PIS";
+        case {"antarctic", "antarctic ice sheet", "ais"}
+            code = "AIS";
+        case {"west antarctic", "west antarctic ice sheet", "wais"}
+            code = "WAIS";
+        case {"east antarctic", "east antarctic ice sheet", "eais"}
+            code = "EAIS";
+        otherwise
+            error("Unrecognised ice sheet name: %s", name)
+    end
+
+end
+
+function isIcesheet = labelIceSheet(lonlat, icesheet)
+
+    arguments (Input)
+        lonlat (:, 2) {mustBeNumeric}
+        icesheet {mustBeTextScalar}
+    end
+
+    arguments (Output)
+        isIcesheet (:, 1) logical
+    end
+
+    isIcesheet = isinterior(icesheetPoly(icesheet), lonlat(:, 1), lonlat(:, 2));
+
+end