Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/README.md =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/README.md (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/README.md (revision 26975) @@ -0,0 +1,15 @@ +# Function List + +* analyzeCalvingFront.m -- Interpolate the transient solutions to the calving front positions of the given flowlines +* averageOverTime.m -- For a given range, average transient data in time +* computeGrad.m -- Compute the gradient of a given field +* extractTransientSolutions.m -- Load a model with transient solutions, and put each of the items into individual array +* interpZeroPos.m -- Find zero positions (y=0) of the given curve defined by (X,Y) +* integrateOverDomain.m -- integrate a given variable over the whole domain, can set a mask +* integrateOverTime.m -- integrate a time series over the given time steps +* inverseHistograms.m -- Compute the inverse quantile function +* npsd.m -- Compute normalized power spectral density +* projectToFlowlines.m -- Project data from a mesh to a given flowline +* psd.m -- Compute power spectral density +* rescalegradientNan.m -- Replace Nan by 0 in the cost function, then rescale it +* wfDistance.m -- Compute Wasserstein-Fourier distance Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/analyzeCalvingFront.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/analyzeCalvingFront.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/analyzeCalvingFront.m (revision 26975) @@ -0,0 +1,62 @@ +function [flowline, icemaskFL] = analyzeCalvingFront(md, flowline, transientSolutions) +%analyzeCalvingFront - to find calving front position at the given flowline, +% if the given model contains transient solutions, then +% return a time dependent calving front with +% [Xmain, positionx, positiony, time], icemask, calving rate, +% melting rate and velocity magnitude, sigma, thickness +% Author: Cheng Gong +% Last modified: 2020-09-25 + +transient = isfield(md.results,'TransientSolution'); + +if transient + % extract data from model + icemaskFL = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,transientSolutions.ice_levelset,flowline.x,flowline.y); + + % solve for the calving front coordinates + cf = interpZeroPos([flowline.Xmain(:), flowline.x, flowline.y], icemaskFL); + positionx = cf(:, 2); + positiony = cf(:, 3); + + % Calving rate C + cRate = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,... + transientSolutions.calvingRate,positionx, positiony); + cRate = diag(cRate); % not very efficient, but works + + % melting rate M + mRateTemp = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,... + transientSolutions.meltingRate,positionx, positiony); + mRate = diag(mRateTemp); % not very efficient, but works + + % velocity + vel = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,... + transientSolutions.vel, positionx, positiony); + vel = diag(vel); + + % sigmaVM + SigmaVM = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,... + transientSolutions.SigmaVM, positionx, positiony); + SigmaVM = diag(SigmaVM); + + % calving front + calvingFront = [cf, transientSolutions.time(:)]; + + % solutions along the flowline + SigmaVM_FL = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,... + transientSolutions.SigmaVM,flowline.x, flowline.y); + + % Thickness along the flowlines + thickness = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,... + transientSolutions.thickness, flowline.x, flowline.y); + + % save to flowline + flowline.calvingFront = calvingFront; + flowline.cRate = cRate; + flowline.mRate = mRate; + flowline.vel = vel; + flowline.SigmaVM = SigmaVM; + flowline.SigmaVM_FL = SigmaVM_FL; + flowline.thickness = thickness; +else + error('Calving front detectiong for the steady state is not implemented yet!') +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/averageOverTime.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/averageOverTime.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/averageOverTime.m (revision 26975) @@ -0,0 +1,86 @@ +function averagedData = averageOverTime(data, time, startP, endP) +%averageOverTime - compute the time averaged value of data in the range [startP, endP] +% +% data: time dependent data, +% time: has same number of columns as data, sorted in ascending order +% startP: start time point for averaging, can be at any time point, or +% even not coincide with the given points in 'time' +% endP: end time point for averaging. if this is the same as the startP, +% or not given, then the output will be the linear interpolation of the +% data at startP. +% +% if out of the given time range, use constant extrapolation, +% if within the range, do linear interpolation. +% +% Author: Cheng Gong +% Last modified: 2020-09-09 +if nargin < 4 + endP = startP; +end + +Nt = length(time); +if Nt ~= size(data,2) + error('The data and time does not match!'); +end + +if startP > endP + error('The given start point has to be smaller than the end point!'); +elseif startP == endP + % extrapolation if earlier than the first time point + if startP <= time(1) + i1 = 1; + i2 = 2; + % extrapolation if later than the last time point + elseif startP >= time(Nt) + i1 = Nt-1; + i2 = Nt; + else + pos = (time > startP); + i2 = find(pos, 1, 'first'); + i1 = i2 - 1; + end + % linear interpolation + averagedData = data(:, i1) + (data(:, i1)-data(:, i2)).*(startP-time(i1))./(time(i1)-time(i2)); +else + % find the first and last ID in the time series within the given range + pos = ((time>=startP) & (time<=endP)); + firstId = find(pos, 1, 'first'); + lastId = find(pos, 1, 'last'); + if isempty(firstId) | isempty(lastId) + averagedData = nan(size(data,1), 1); + % error('Time points out of range'); + else + % computed the integral with trapzoidal rule + integData = zeros(size(data,1), 1); + for i = firstId:lastId-1 + integData = integData + (data(:,i+1) + data(:,i)) .* (time(i+1) - time(i)).*0.5; + end + + % special treatment for the first and last time step, which are not + % complete steps. + if firstId == 1 + % extrapolation + integData = integData + data(:,1) .* (time(1)-startP); + else + % interpolation + integData = integData + 0.5.*(time(firstId)-startP)./ ... + (time(firstId)-time(firstId-1)) .* ... + ((time(firstId)-startP) .* data(:,firstId-1) + ... + (time(firstId)+startP-2*time(firstId-1)).*data(:,firstId)); + end + + if lastId == Nt + % extrapolation + integData = integData + data(:,Nt) .* (endP-time(Nt)); + else + % interpolation + integData = integData + 0.5.*(time(lastId)-endP)./ ... + (time(lastId+1)-time(lastId)) .* ... + ((time(lastId)-endP) .* data(:,lastId+1) + ... + (time(lastId)+endP-2*time(lastId+1)).*data(:,lastId)); + end + + averagedData = integData ./ (endP-startP); + end +end + Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m (revision 26975) @@ -0,0 +1,45 @@ +function [gradx, grady]=computeGrad(index,x,y,field) +%COMPUTEHESSIAN - compute the gradient from a field + +%some variables +numberofnodes=length(x); +numberofelements=size(index,1); +numberoftime = size(field, 2); + +%some checks +if (length(field)~=numberofnodes) && (length(field)~=numberofelements) + error('ComputeHessian error message: the given field size not supported yet'); +end +%initialization +line=index(:); +linesize=3*numberofelements; + +%get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma +[alpha, beta]=GetNodalFunctionsCoeff(index,x,y); +areas=GetAreas(index,x,y); + +%compute weights that hold the volume of all the element holding the node i +weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofnodes,1); + +%compute field on nodes if on elements +if length(field)==numberofelements + field=sparse(line,ones(linesize,1),repmat(areas.*field,3,1),numberofnodes,1)./weights ; +end + +%Compute gradient for each element +if numberoftime == 1 + grad_elx=sum(field(index).*alpha,2); + grad_ely=sum(field(index).*beta,2); +else + grad_elx = zeros(numberofelements,numberoftime); + grad_ely = zeros(numberofelements,numberoftime); + for i = 1:3 + grad_elx = grad_elx + field(index(:,i), :).*alpha(:,i); + grad_ely = grad_ely + field(index(:,i), :).*beta(:,i); + end +end +%Compute gradient for each node (average of the elements around) +gradx=sparse(repmat(line,1,numberoftime),cumsum(ones(linesize,numberoftime),2),repmat(areas.*grad_elx,3,1),numberofnodes,numberoftime); +grady=sparse(repmat(line,1,numberoftime),cumsum(ones(linesize,numberoftime),2),repmat(areas.*grad_ely,3,1),numberofnodes,numberoftime); +gradx=gradx./weights; +grady=grady./weights; Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m (revision 26975) @@ -0,0 +1,14 @@ +function [transientSolutions]=extractTransientSolutions(md) +%extractTransientSolutions - take the transient solutions out from model +% and put each of them into an individual array. +transientSolutions.time = cell2mat({md.results.TransientSolution(:).time}); +transientSolutions.vx = cell2mat({md.results.TransientSolution(:).Vx}); +transientSolutions.vy = cell2mat({md.results.TransientSolution(:).Vy}); +transientSolutions.vel = cell2mat({md.results.TransientSolution(:).Vel}); +transientSolutions.volume = cell2mat({md.results.TransientSolution(:).IceVolume}); +transientSolutions.thickness = cell2mat({md.results.TransientSolution(:).Thickness}); +transientSolutions.SigmaVM = cell2mat({md.results.TransientSolution(:).SigmaVM}); +transientSolutions.smb = cell2mat({md.results.TransientSolution(:).SmbMassBalance}); +transientSolutions.ice_levelset = cell2mat({md.results.TransientSolution(:).MaskIceLevelset}); +transientSolutions.calvingRate = cell2mat({md.results.TransientSolution(:).CalvingCalvingrate}); +transientSolutions.meltingRate = cell2mat({md.results.TransientSolution(:).CalvingMeltingrate}); Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/fillInNan.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/fillInNan.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/fillInNan.m (revision 26975) @@ -0,0 +1,26 @@ +function newdata = fillInNan(md, data, varargin) +%fillInNan - use the mean of surrouding data to iteratidatay fill in the Nan in data +% data should have the same size as md.mesh.x +% +% Author: Cheng Gong +% Last modified: 2021-12-08 +nanvflag = find(isnan(data)); +NNanv = length(nanvflag); + +options = pairoptions(varargin{:}); +maxiter = getfieldvalue(options,'maxiter', 10); +count = 1; + +while((NNanv>0) & (count<=maxiter)) + disp(['Iter ', num2str(count), ': found ', num2str(NNanv), ' Nan in the initial data, fill in them with the mean of their surroudings.']); + for i = 1:NNanv + [eleid, ~] = find(md.mesh.elements == nanvflag(i)); + nodeid = md.mesh.elements(eleid, :); + data(nanvflag(i)) = mean(data(nodeid(:)) , 'omitnan'); + end + nanvflag = find(isnan(data)); + NNanv = length(nanvflag); + count = count + 1; +end +newdata = data; + Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m (revision 26975) @@ -0,0 +1,32 @@ +function [intData, meanData, areas] = integrateOverDomain(md, data, masked, weights) +% integrateOverDomain - integrating data over the whole domain +% +% intData: integral of the data over each element +% meanData: intData/areas +% areas: areas of the domain +if nargin < 4 + weights = ones(size(data)); + if nargin<3 + masked = []; + end +end + +% Set the area with masked=1 to nan +data(masked) = nan; +weights(masked) =nan; + +% get the mesh +elements=md.mesh.elements; +x=md.mesh.x; +y=md.mesh.y; + +%compute areas; +eleAreas=GetAreas(elements,x,y); + +% integrate nodal data to element +eleData = 1/3*eleAreas.*(data(elements(:,1),:).*weights(elements(:,1),:) + data(elements(:,2),:).*weights(elements(:,2),:) + data(elements(:,3),:).*weights(elements(:,3),:)); +eleAreas = 1/3*eleAreas.*(weights(elements(:,1),:)+weights(elements(:,2),:)+weights(elements(:,3),:)); + +intData = sum(eleData(:),'omitnan'); +areas = sum(eleAreas(:),'omitnan'); +meanData = intData / areas; Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverTime.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverTime.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverTime.m (revision 26975) @@ -0,0 +1,26 @@ +function [newData, intData] = integrateOverTime(tdata, data, time, Nintdt) +% integrateOverDomain - integrating data over the whole domain +% +% Input: +% tdata: time steps of data +% data: time series of the data +% time: time steps to be used for integration +% Nintdt: number of time steps in 'time' to be used for integration +% +% Return: +% newData: data projected on to 'time' +% intData: integrated 'newData' + +dt = abs(time(2) -time(1)); +% patch the first N time steps +intTime = [[-Nintdt:-1]*dt+time(1), time]; +% integration scheme +integ = ones(1, Nintdt) * dt; + +% interpolate data from the time scheme 'tdata' to 'time' +newData = interp1(double(tdata), double(data), intTime, 'linear', 'extrap'); + +% integrated, only works for equidistance time steps for now +intData = conv(integ, newData); +intData = intData(Nintdt:end-Nintdt); +newData = newData(Nintdt+1: end); Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/interpZeroPos.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/interpZeroPos.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/interpZeroPos.m (revision 26975) @@ -0,0 +1,52 @@ +function X0 = interpZeroPos(X, Y) +%interpZeroPos - find the x values for y=0 with line interpolation of the +% given X and Y coordinates. +% X: is a matrix of size mxn, each column represent a dimension +% Y: is a matrix of size mxk, k is the number of date sets, and each +% data set must contains positive and negative values. +% The return value: +% x0: is a matrix of size kxn, the rows correspond to each column +% of X and the columns are for each data set. +% +% The method is to find the two consective values y1 and y2 where +% Y changes the sign, and the corresponding X are x1 and x2. +% x2*y1-x1*y2 +% x0 = ------------- +% y1-y2 +% Author: Cheng Gong +% Last modified: 2020-10-20 +[m, n] = size(X); +[my, k] = size(Y); +if m~=my + error('rows of X and Y are not the same!'); +end + +% find the place where Y changes signs, or y=0 at some points +mask = (Y(1:m-1,:).*Y(2:m,:) <=0); +[row, col] = find(mask); +ind = row; + +% use only the first index from each column +if length(col)>k + disp('Multiple candidates of the y=0 values are found in one data set! Be defalut, take the first y=0 from upstream.'); + [col, ia, ~] = unique(col); + ind = row(ia); +end + +if length(col)epsilon); +Sa = S(A); + +% cummulative sum +cdf = cumtrapz(S, mu); +cdfa = cdf(A); + +% set the first value to 0 and last value to 1 +cdfa(1) = 0; +Sa(1) = 0; + +if cdfa(end) < 1 + cdfa = [cdfa(:);1]; + Sa = [Sa(:);S(end)]; +end + +[~, ind] = unique(cdfa); +% linear interpolation +q_Sinv = interp1(cdfa(ind), Sa(ind), Sinv); + +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/movingAverage.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/movingAverage.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/movingAverage.m (revision 26975) @@ -0,0 +1,42 @@ +function outdata = movingAverage(rawdata, varargin) +% +% rawdata is an MxN matrix, where the first M-1 rows are the raw data, the row M is the time + +options = pairoptions(varargin{:}); +timeWindow = getfieldvalue(options, 'time window', 30); % unit in days +resamp = getfieldvalue(options, 'resample', 1); % resample the data set with the size of window + +% get data info {{{ +[M, N] = size(rawdata); +if M < 2 + error('The data has to be at least two rows!'); +end +if N < 2 + error('The data has to be at least two columns!'); +end + +% sort the time +time = rawdata(M, :); +[time, tI] = sort(time); +data = rawdata(1:M-1, :); +data = data(:, tI); + +dt = time(2) - time(1); +if dt < 0 + error('The time need to be in ascending order.'); +end +windowSize = ceil(timeWindow/365/dt); +%}}} +%% filter {{{ +if resamp + ind = [1:timeWindow:N]; +else + ind = [1:N]; +end +b = (1/windowSize)*ones(1,windowSize); +a = 1; +disp(['Applying filter of width ', num2str(timeWindow)]); +smoothdata = filter(b,a,data,[],2); +time = time(ind); +outdata = [smoothdata(:,ind);time]; +%}}} Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/npsd.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/npsd.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/npsd.m (revision 26975) @@ -0,0 +1,32 @@ +function [npsdx, psdx, freq] = npsd(x) +%npsd - Normalized Power Spectral Density (using fft) +% +% x: time series +% +% psdx: Power Spectral Density(right half plane, value doubled) +% freq: frequency(non-negative half) +% +% Author: Cheng Gong +% Date: 2021-08-17 + +% length of data +Nx = length(x); +% fft +xdft = fft(x); +xdft = xdft(1:floor(Nx/2+1)); +% power specturm +psdx = (1/(2*pi*Nx)) * abs(xdft).^2; +% In order to conserve the total power, multiply all frequencies that +% occur in both sets - the positive and negative frequencies — by a factor of 2 +psdx(2:end-1) = 2*psdx(2:end-1); +freq = 0:(2*pi)/Nx:pi; +% normalize(trapzoidal rule) +intS = trapz(freq, psdx); +if intS == 0 + psdx = 1 +psdx; + intS = trapz(freq, psdx); +end +npsdx = psdx ./ intS; + + +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/projectToFlowlines.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/projectToFlowlines.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/projectToFlowlines.m (revision 26975) @@ -0,0 +1,16 @@ +% Function to project values on the mesh to a flowline +% +% md - ISSM model with mesh +% pValue - data on the mesh +% fx - x coordinates of the flowline +% fy - y coordinates of the flowline +% +% Author: Cheng Gong +% Last modified: 2021-01-27 + +function valueC = projectToFlowlines(md, pValue, fx, fy) + temp = InterpFromMeshToMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y,... + pValue, fx, fy); + temp = diag(temp); % not very efficient, but works + valueC = temp(:)'; +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/psd.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/psd.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/psd.m (revision 26975) @@ -0,0 +1,24 @@ +function [psdx, freq] = psd(x) +%psd - Power Spectral Density (using fft) +% +% x: time series +% +% psdx: Power Spectral Density(right half plane, value doubled) +% freq: frequency(non-negative half) +% +% Author: Cheng Gong +% Date: 2021-08-17 + +% length of data +Nx = length(x); +% fft +xdft = fft(x); +xdft = xdft(1:floor(Nx/2+1)); +% power specturm +psdx = (1/(2*pi*Nx)) * abs(xdft).^2; +% In order to conserve the total power, multiply all frequencies that +% occur in both sets - the positive and negative frequencies — by a factor of 2 +psdx(2:end-1) = 2*psdx(2:end-1); +freq = 0:(2*pi)/Nx:pi; + +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/rescalegradientNan.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/rescalegradientNan.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/rescalegradientNan.m (revision 26975) @@ -0,0 +1,7 @@ +% function to replace Nan by 0 in J, then rescale with the mass matrix +% +function J = rescalegradientNan(md, unscaledJ) + +nanflag = isnan(unscaledJ); +unscaledJ(nanflag) = 0; +J = rescalegradient(md, unscaledJ); Index: /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/wfDistance.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/wfDistance.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/wfDistance.m (revision 26975) @@ -0,0 +1,26 @@ +function w2 = wfDistance(x, y) +%wfDistance - compute Wasserstein-Fourier Distance, according to Cazelles, E., +% Robert, A. & Tobar, F. The Wasserstein-Fourier Distance for +% Stationary Time Series. Ieee T Signal Proces 69, 709–721 (2019). +% +% x: first time series +% y: second time series +% +% w2: Wasserstein-Fourier Distance, or W2 distance of sx and sy +% +% Author: Cheng Gong +% Date: 2021-08-18 +N = length(x); +Sinv=linspace(0,1,N); + +[csx, ~, freqx] = npsd(x); +[csy, ~, freqy] = npsd(y); + +% inverse cummulative function +[~, qx] = inverseHistograms(csx, freqx, Sinv); +[~, qy] = inverseHistograms(csy, freqy, Sinv); + +% +w2 = sqrt(trapz(Sinv ,(qx-qy).^2)); + +end Index: /issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m (revision 26975) @@ -0,0 +1,69 @@ +function expxy2shpll(expfilename,shpfilename,geometry) +%expxy2shpll +% Convert .exp to .shp file +% +% Usage: +% expxy2shpll('glacier.exp', 'glacier.shp') +% expxy2shpll('glacier.exp', 'glacier.shp', geometry) +% +% geometry (optional)- 'MultiPoint' : point clouds +% 'Point' : single point +% 'Line' : two points line +% 'polygon' : multiple points +% + +%check file extensions +[pathstr,name,ext] = fileparts(shpfilename); +if ~strcmp(ext,'.shp'), + error(['Shapefile ' shpfilename ' does not have an extension .shp']); +end + +[pathstr,name,ext] = fileparts(expfilename); +if ~strcmp(ext,'.exp'), + error(['Exp file ' expfilename ' does not have an extension .exp']); +end + +shp=expread(expfilename); +[lat, lon] = xy2ll(shp.x, shp.y, 1, 45, 70); +shp.x = lon; +shp.y = lat; + +%initialize number of profile +count=1; + +contours=struct([]); +for i=1:length(shp), + if nargin < 3 + + %TEMP + %if contains(shp(i).name,'_pointcloud'); + % continue; + %end + + if length(shp(i).x) == 0 + continue; + elseif contains(shp(i).name,'_pointcloud'); + geometry = 'MultiPoint'; + shp(i).name = erase(shp(i).name,'_pointcloud'); + elseif length(shp(i).x) == 1 + geometry = 'Point'; + elseif length(shp(i).x) < 3 + geometry = 'Line'; + else + if (shp(i).x(end)==shp(i).x(1) && shp(i).y(end)==shp(i).y(1)), + geometry = 'Polygon'; + else + geometry = 'Line'; + end + end + end + contours(count).Geometry=geometry; + contours(count).id=i; + contours(count).Name=shp(i).name; + contours(count).X=shp(i).x; + contours(count).Y=shp(i).y; + count = count+1; +end + +%Make sure it is one single geometry otherwise it will yell at you +shapewrite(contours,shpfilename); Index: /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m (revision 26974) +++ /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m (revision 26975) @@ -1,9 +1,33 @@ -% This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable -% tif data in /totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/ within the given time period (in decimal years) -% For some reason, each .tif file in this folder contains two sets of data, only the first dataset is useful +function dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend,varargin) +%interpFromMEaSUREsGeotiff: +% This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable +% tif data in /totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/ within the given time period (in decimal years) +% For some reason, each .tif file in this folder contains two sets of data, only the first dataset is useful +% +% Usage: +% dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend, varargin) +% +% X, Y are the coordinates of the mesh +% Tstart and Tend decimal year of the start and end time +% +% Example: +% obsData = interpFromMEaSUREsGeotiff(md.mesh.x,md.mesh.y, tstart, tend); +% +% Options: +% - 'glacier': which glacier to look for +options = pairoptions(varargin{:}); +glacier = getfieldvalue(options,'glacier','Jakobshavn'); -function dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend) - -foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/'; +if strcmp(glacier, 'Jakobshavn') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/'; +elseif strcmp(glacier, 'Kangerlussuaq') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Kangerlussuaq_2006_2021/'; +elseif strcmp(glacier, 'Store') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Store_2008_2021/'; +elseif strcmp(glacier, 'Rink') + foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Rink_2008_2022/'; +else + error(['The velocity data for ', glacier, ' is not available, please download from NSIDC first.']); +end % get the time info from file names Index: /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpISMIP6Temp.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpISMIP6Temp.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpISMIP6Temp.m (revision 26975) @@ -0,0 +1,15 @@ +function temp = interpISMIP6Temp(X, Y) +%interpISMIP6Temp - interpolate land ice basal temperature to the given mesh +% +% X and Y are the coordinates of the mesh +% +% Author: Cheng Gong +% Last modified: 2021-12-06 + +filename = '/totten_1/ModelData/Greenland/ISMIP6/GreenlandISMIP6-Morlighem-2020-10-01.nc'; + +x = ncread(filename, 'x'); +y = ncread(filename, 'y'); +tb = ncread(filename, 'tb'); + +temp = InterpFromGrid(double(x), double(y), double(tb'), X, Y); Index: /issm/trunk-jpl/src/m/contrib/chenggong/visualization/imageNonUni.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/visualization/imageNonUni.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/visualization/imageNonUni.m (revision 26975) @@ -0,0 +1,23 @@ +function imageNonUni(varargin) +%imageNonUni - draw an image with non uniform grid +% +% Author: Cheng Gong +% Last modified: 2020-08-24 +x = varargin{1}; +y = varargin{2}; +C = varargin{3}; + +nx = length(x); +ny = length(y); +[ncy, ncx] = size(C); + +if ((nx ~= ncx) || (ny~=ncy)) + error('The size of x-y coordinates do not match the data!'); +end + +X = repmat(x(:)',ny , 1); +Y = repmat(y(:), 1, nx); + + +surf(X, Y, zeros(ny,nx), C, 'EdgeColor','None'); +view(2); Index: /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m (revision 26975) +++ /issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m (revision 26975) @@ -0,0 +1,48 @@ +function plotCompareTransientFlowline(velList, titleList, md, flowline, TStart, TEnd, xl, yl, ca) +%plotCompareTransientFlowline - to compare the solutions along a given +% flowine for the given time periods +% +% velList: a cell of transient velocities on the model mesh +% titleList: names of the velocity data +% md: ISSM model +% flowline: a flowline data structure contains flowline.x and flowline.y. +% If possible, it should also contains flowline.Xmain, which is +% the distance along the flowline from upstream to the terminus position. +% TStart: a list of start time points +% TEnd: a list of end time points +% ca: caixs value in plot +% xl: xlim value +% yl: ylim value +% +% Author: Cheng Gong +% Date: 2021-12-06 + +N = length(velList); +vel_flowline = cell(N,1); +if ~isfield(flowline, 'Xmain') + flowline.Xmain = cumsum([0;sqrt((flowline.x(2:end)- flowline.x(1:end-1)).^2 + (flowline.y(2:end)- flowline.y(1:end-1)).^2)]')/1000; +end + +for i = 1:N + vel_flowline{i} = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,velList{i},flowline.x,flowline.y); +end + +for p = 1:N + subplot(N,1,p) + hold on + for i = 1:length(TStart) + overlay = [vel_flowline{p}(:,i), vel_flowline{p}(:,i)]'; + imAlpha = ones(size(overlay)); + imAlpha(isnan(overlay)) = 0; + imageNonUni(flowline.Xmain, [TStart(i), TEnd(i)], overlay); + end + title(titleList{p}, 'Interpreter', 'latex'); + xlim(xl) + xlabel('x(km)') + ylabel('year') + ylim(yl); + caxis(ca) + colormap('jet') + h=colorbar; + title(h,'m/a') +end