Skip to content

Commit e143be9

Browse files
hvandermarelhvandermarel
committed
added itrfmap, with demo and live editor demo.
1 parent 6efee3f commit e143be9

File tree

3 files changed

+390
-0
lines changed

3 files changed

+390
-0
lines changed

‎demoitrfmap.m‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
%% Plot map of ITRF and ETRF velocity or coordinate changes
2+
%
3+
% Script to plot ITRF and ETRF velocity and position differences, with
4+
% or without plate motion model correction.
5+
%
6+
% Created: 9 February 2018 by Hans van der Marel
7+
% Modified: 2 june 2025 by Hans van der Marel
8+
% - Major overhaul using itrf toolbox version 2 (with newly
9+
% added pmm functions and support for more TRF's)
10+
% - optional correction for plate motion (pmm)
11+
% - choice between plotting velocity differences and position
12+
% differences
13+
% - added option for low-resolution maps (useful when high
14+
% resolution shape files are unavailable)
15+
16+
%% Add required toolboxes to the Matlab path
17+
%
18+
% Requires the Matlab mapping toolbox (if not available, should not be
19+
% too difficult to adjust the plotting part using coastlines only).
20+
21+
addtoolbox('itrf');
22+
addtoolbox('crsutil'); % import plh2xyz and xyz2neu from crsutil
23+
24+
% addpath('d:\Surfdrive\Matlab\toolbox\itrf')
25+
% addpath('d:\Surfdrive\Matlab\toolbox\crsutil')
26+
27+
%% Legend to itrfvelmap arguments
28+
%
29+
% Reference frames Plate motion model (none/pmm) -----+
30+
% from ------------+ |
31+
% to ------+ | Plot type (vel/pos) --+ |
32+
% | | | |
33+
% v v v v
34+
% itrfmap(refto, reffrom, refepoch, latrange, lonrange, sel, pmm, tecplate, id)
35+
36+
%% Velocity and position differences over Europe
37+
38+
% ROI and map limits for Europe
39+
40+
latrange=[35 68]; lonrange=[-12 35 ];
41+
42+
% Examples of itrfvelmap (uncomment one or more)
43+
44+
% itrfmap('ITRF2008', 'ETRF2000', 2010.0, latrange, lonrange)
45+
% itrfmap('ITRF2014', 'ETRF2000', 2018.0, latrange, lonrange, 'vel')
46+
% itrfmap('ITRF2020', 'ETRF2000', 2020.0, latrange, lonrange, 'vel', 'none')
47+
itrfmap('ITRF2008', 'ETRF2000', 2020.0, latrange, lonrange, 'vel', 'GSRM 2.1 IGS08', 'Eurasia')
48+
%
49+
% itrfmap('ITRF2020', 'ITRF2014', 2020.0, latrange, lonrange, 'vel', 'none')
50+
%
51+
% itrfmap('ITRF2008', 'ETRF2020', 1989.0, latrange, lonrange, 'pos', 'none')
52+
% itrfmap('ITRF2008', 'ETRF2020', 2020.0, latrange, lonrange, 'pos', 'GSRM 2.1 IGS08', 'Eurasia')
53+
% itrfmap('ITRF2020', 'ETRF2020', 1989.0, latrange, lonrange, 'pos', 'none')
54+
% itrfmap('ITRF2020', 'ETRF2020', 2020.0, latrange, lonrange, 'pos', 'ITRF2020', 'Eurasia')
55+
56+
57+
%% Velocity and position differences over the Netherlands
58+
59+
% Add optional folder with high-resolution land areas, rivers and borders
60+
% to Matlab path, if not available, low resolution maps that come with the
61+
% Matlab mapping toolbox are used.
62+
63+
addpath('d:\Surfdrive\Matlab\models\')
64+
65+
% ROI and map limits for the Netherlands
66+
67+
latrange=[50 54.5];
68+
lonrange=[1.5 8.5 ];
69+
70+
% Examples of itrfvelmap (uncomment one or more)
71+
72+
itrfmap('ITRF2008', 'ETRF2000', 2020, latrange, lonrange, 'vel', 'GSRM 2.1 IGS08', 'Eurasia','nweurope')
73+
% itrfmap('ETRF2020', 'ETRF2014', 2020, latrange, lonrange, 'vel', 'none', 'none', 'nweurope')
74+

‎itrfmap.m‎

Lines changed: 316 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,316 @@
1+
function itrfmap(refto, reffrom, refepoch, latrange, lonrange, sel, pmm, tecplate,id)
2+
%ITRFMAP Plot map of ITRF/ETRF velocity and coordinate changes.
3+
% ITRFMAP(REFTO,REFFROM,EPOCH,LATRANGE,LONGRANGE) plots the difference
4+
% in velocity for the reference frame REFTO with the velocity in REFFROM
5+
% on a map with latitude range LATRANGE and longitude range LONRANGE.
6+
% EPOCH is the epoch, but this is only used internally for the computation
7+
% and has no significance for the velocity output.
8+
%
9+
% ITRFMAP(...,SEL), with SEL='pos', plots the coordinate differences
10+
% instead of the velocity differences (SEL='vel' is the default). EPOCH
11+
% is the epoch for the coordinate difference.
12+
%
13+
% ITRFMAP(...,SEL,PMM,TECPLATE) applies a correction to REFTO using
14+
% the plate motion model given by PMM for tectonic plates TECPLATE.
15+
%
16+
% ITRFMAP(...,SEL,PMM,TECPLATE,ID) reads high-resolution shape files
17+
% for land areas, rivers and borders (if available) with ID. The only
18+
% high-resolution data currently supported is ID='nweurope` and the
19+
% shapefiles must on the search path. Default (ID='none') is to use
20+
% low resolution maps that come with the Matlab mapping toolbox.
21+
%
22+
% A short legend to the available arguments is:
23+
%
24+
% Reference frames Plate motion model (none/pmm) -----+
25+
% from ------------+ |
26+
% to ------+ | Plot type (vel/pos) --+ |
27+
% | | | |
28+
% v v v v
29+
% itrfmap(refto, reffrom, epoch, latrange, lonrange, sel, pmm, tecplate, id)
30+
%
31+
% Supported reference frames are [ "ETRF2020"; "ETRF2014"; "ETRF2000";
32+
% "ITRF2020"; "ITRF2014"; "ITRF2008"; "ITRF2005"; "ITRF2000"; "ITRF97"; "ITRF96";
33+
% "ITRF94"; "ITRF93"; "ITRF92"; "ITRF91"; "ITRF90"; "ITRF89"; "ITRF88" ].
34+
%
35+
% Supported plate motion models PMM are [ "none"; "NUVEL 1A NNR"; "ITRF2000";
36+
% "ITRF2008"; "ITRF2014"; "ITRF2020"; "GSRM 2.1 IGS08"; "GSRM 2.1 NNR" ].
37+
% The tectonic plate TECPLATE can be specified by it's full name, four-
38+
% and two character abbreviation. The supported plates depend on the
39+
% plate motion model being used. Some of the possible values for TECPLATE
40+
% are [ "Africa"; "Antarctica"; "Arabia"; "Australia"; "Caribbean"; "Cocos";
41+
% "Eurasia"; "India"; "Nazca"; "NorthAmerica"; "Pacific"; "SouthAmerica" ];
42+
%
43+
% Examples:
44+
%
45+
% latrange=[35 68]; lonrange=[-12 35 ]; % Europe
46+
%
47+
% itrfmap('ITRF2008', 'ETRF2000', 2010.0, latrange, lonrange)
48+
% itrfmap('ITRF2014', 'ETRF2000', 2018.0, latrange, lonrange, 'vel')
49+
% itrfmap('ITRF2020', 'ETRF2000', 2020.0, latrange, lonrange, 'vel', 'none')
50+
% itrfmap('ITRF2008', 'ETRF2000', 2020.0, latrange, lonrange, 'vel', 'GSRM 2.1 IGS08', 'Eurasia')
51+
%
52+
% itrfmap('ITRF2020', 'ITRF2014', 2020.0, latrange, lonrange, 'vel', 'none')
53+
%
54+
% itrfmap('ITRF2008', 'ETRF2020', 1989.0, latrange, lonrange, 'pos', 'none')
55+
% itrfmap('ITRF2008', 'ETRF2020', 2020.0, latrange, lonrange, 'pos', 'GSRM 2.1 IGS08', 'Eurasia')
56+
% itrfmap('ITRF2020', 'ETRF2020', 1989.0, latrange, lonrange, 'pos', 'none')
57+
% itrfmap('ITRF2020', 'ETRF2020', 2020.0, latrange, lonrange, 'pos', 'ITRF2020', 'Eurasia')
58+
%
59+
% latrange=[50 54.5]; lonrange=[1.5 8.5 ]; % Netherlands
60+
% addpath('d:\Surfdrive\Matlab\models\') % Path to high resolution maps
61+
%
62+
% itrfmap('ITRF2008', 'ETRF2000', 2020.0, latrange, lonrange, 'vel', 'GSRM 2.1 IGS08', 'Eurasia','nweurope')
63+
% itrfmap('ETRF2020', 'ETRF2014', 2020.0, latrange, lonrange, 'vel', 'none', 'none', 'nweurope')
64+
%
65+
% See also ITRF2ITRF and PMMVEL.
66+
%
67+
% Requires the Matlab mapping toolbox (if not available, should not be
68+
% too difficult to adjust the plotting part using coastlines only).
69+
%
70+
% (c) Hans van der Marel, Delft University of Technology, 2018-2025.
71+
72+
% Created: 9 February 2018 by Hans van der Marel
73+
% Modified: 2 June 2025 by Hans van der Marel
74+
% - Major overhaul using itrf toolbox version 2 (with newly
75+
% added pmm functions and support for more TRF's)
76+
% - optional correction for plate motion (pmm)
77+
% - choice between plotting velocity differences and position
78+
% differences
79+
% - added option for low-resolution maps (useful when high
80+
% resolution shape files are unavailable)
81+
% 14 July 2025 by Hans van der Marel
82+
% - single figure with subplots
83+
% - added extra argument for high resolution maps
84+
% - added function help
85+
86+
tic
87+
88+
% Check arguments
89+
90+
if nargin < 9
91+
id='none';
92+
if nargin < 7
93+
pmm='none';
94+
if nargin < 6
95+
sel='vel';
96+
end
97+
end
98+
end
99+
if nargin < 5 || ( nargin == 7 && ~strcmpi(pmm,'none') )
100+
error('insufficient number of input arguments')
101+
end
102+
103+
% Compute meshgrid
104+
105+
bbox = [ lonrange' latrange'];
106+
107+
%opt.interval=.1; % interpolation interval in deg
108+
opt.interval= (ceil(lonrange(2))-floor(lonrange(1))+ceil(latrange(2))-floor(latrange(1)))/100;
109+
tx=floor(bbox(1)):opt.interval:ceil(bbox(2));
110+
ty=floor(bbox(3)):opt.interval:ceil(bbox(4));
111+
112+
[qx,qy]=meshgrid(tx,ty);
113+
114+
% Compute (nominal) velocity/position REFTO wrt to REFFROM
115+
116+
plh=[ qy(:)*pi/180 qx(:)*pi/180 zeros(size(qx(:)))];
117+
xyz=plh2xyz(plh);
118+
119+
xyzitrf=itrf2itrf([xyz zeros(size(xyz)) ] ,char(reffrom),char(refto),refepoch);
120+
121+
posneuitrf=xyz2neu(xyzitrf(:,1:3)-xyz,xyz,'r');
122+
velneuitrf=xyz2neu(xyzitrf(:,4:6),xyz,'r');
123+
124+
% Optionally correct for plate motion model (PMM) velocity
125+
126+
if ~strcmpi(pmm,'none')
127+
posneuitrf = posneuitrf - pmmvel(pmm,tecplate,plh) * ( refepoch - 1989.0); %ETRS89 ref epoch
128+
velneuitrf = velneuitrf - pmmvel(pmm,tecplate,plh);
129+
end
130+
131+
% Select plot type
132+
133+
if strncmpi(sel,'velocity',3)
134+
plttitle = [ 'Velocity{\Delta} ' char(refto) ' wrt ' char(reffrom) ];
135+
units='[mm/y]';
136+
elseif strncmpi(sel,'position',3)
137+
velneuitrf=posneuitrf;
138+
plttitle = [ 'Coordinate{\Delta} ' char(refto) '@' num2str(refepoch,'%.1f') ' wrt ' char(reffrom) ];
139+
units='[mm]';
140+
else
141+
error('Invalid selection, use vel or pos')
142+
end
143+
if ~strcmpi(pmm,'none')
144+
plttitle = [ plttitle ', pmm=' char(pmm)];
145+
end
146+
147+
% Land areas, rivers and borders
148+
149+
if ~strcmpi(id,'none')
150+
try
151+
%id='nweurope';
152+
resolution='f';
153+
landareas = shaperead(['landareas_' id '_' resolution],'UseGeoCoords', true, 'BoundingBox', bbox);
154+
borders = shaperead(['borders_' id '_' resolution],'UseGeoCoords', true, 'BoundingBox', bbox);
155+
rivers = shaperead(['rivers_' id '_' resolution],'UseGeoCoords', true, 'BoundingBox', bbox);
156+
disp('Using high resolution land areas, rivers and borders.')
157+
catch
158+
disp('No high resolution land areas, rivers and borders available, maybe you did not set the path, using low resolution instead.')
159+
id='none';
160+
end
161+
end
162+
if strcmpi(id,'none')
163+
landareas = readgeotable("landareas.shp");
164+
rivers = readgeotable("worldrivers.shp");
165+
borders = geopoint(); % empty
166+
end
167+
168+
% Plot four figures, with basemap and contour overlay of velocities or
169+
% differences, in resp. North, East, Horizontal and Vertical direction
170+
171+
% Fig 1 - Plot the North velocity/difference
172+
173+
qz=reshape(velneuitrf(:,1)*1000,size(qx));
174+
levels=getlevels(qz);
175+
176+
figure('position',[10,10,1000,800]);
177+
t=tiledlayout(2,2,'TileSpacing','Compact','Padding','Compact');
178+
title(t,plttitle)
179+
180+
nexttile;
181+
%figure;
182+
ax=worldmap(latrange,lonrange);
183+
setm(ax,'GLineStyle',':','GColor','k','FFaceColor','none'); % 'FFaceColor','cyan'
184+
%[cs,h]=contourfm(qy,qx,qz,'LevelStep',levelstep,'LineColor',[0.6 0.6 0.6]);
185+
try
186+
[cs,h]=contourfm(qy,qx,qz,levels,'LineColor',[0.6 0.6 0.6]);
187+
ht = clabelm(cs,h);
188+
set(ht,'Color','white','BackgroundColor','none','FontWeight','normal')
189+
catch
190+
disp('Contourfm failure, skip...')
191+
end
192+
c=colorbar;
193+
c.Label.String=units;
194+
geoshow(landareas, 'FaceColor','none','EdgeColor','black');
195+
geoshow(rivers, 'Color', 'cyan')
196+
geoshow(borders, 'Color','black')
197+
%title(['North ' plttitle])
198+
title('North')
199+
200+
% Fig 2 - Plot the East velocity/difference
201+
202+
qz=reshape(velneuitrf(:,2)*1000,size(qx));
203+
levels=getlevels(qz);
204+
205+
nexttile;
206+
%figure;
207+
ax=worldmap(latrange,lonrange);
208+
setm(ax,'GLineStyle',':','GColor','k','FFaceColor','none'); % 'FFaceColor','cyan'
209+
%[cs,h]=contourfm(qy,qx,qz,'LevelStep',levelstep,'LineColor',[0.6 0.6 0.6]);
210+
try
211+
[cs,h]=contourfm(qy,qx,qz,levels,'LineColor',[0.6 0.6 0.6]);
212+
ht = clabelm(cs,h);
213+
set(ht,'Color','white','BackgroundColor','none','FontWeight','normal')
214+
catch
215+
disp('Contourfm failure, skip...')
216+
end
217+
c=colorbar;
218+
c.Label.String=units;
219+
geoshow(landareas, 'FaceColor','none','EdgeColor','black');
220+
geoshow(rivers, 'Color', 'cyan')
221+
geoshow(borders, 'Color','black')
222+
%title(['East ' plttitle])
223+
title('East')
224+
225+
% Fig 3 - Plot the Horizontal velocity/difference
226+
227+
qz=reshape(sqrt(velneuitrf(:,1).^2+velneuitrf(:,2).^2)*1000,size(qx));
228+
levels=getlevels(qz);
229+
230+
nexttile;
231+
%figure;
232+
ax=worldmap(latrange,lonrange);
233+
setm(ax,'GLineStyle',':','GColor','k','FFaceColor','none'); % 'FFaceColor','cyan'
234+
%[cs,h]=contourfm(qy,qx,qz,'LevelStep',levelstep,'LineColor',[0.6 0.6 0.6]);
235+
try
236+
[cs,h]=contourfm(qy,qx,qz,levels,'LineColor',[0.6 0.6 0.6]);
237+
ht = clabelm(cs,h);
238+
set(ht,'Color','white','BackgroundColor','none','FontWeight','normal')
239+
catch
240+
disp('Contourfm failure, skip...')
241+
end
242+
c=colorbar;
243+
c.Label.String=units;
244+
geoshow(landareas, 'FaceColor','none','EdgeColor','black');
245+
geoshow(rivers, 'Color', 'cyan')
246+
geoshow(borders, 'Color','black')
247+
248+
istep=6;
249+
ix=floor(istep/2)+1:istep:size(qx,1);
250+
iy=floor(istep/2)+1:istep:size(qx,2);
251+
qxx=qx(ix,iy);
252+
qyy=qy(ix,iy);
253+
qzv=reshape(velneuitrf(:,1)*1000,size(qx));
254+
qzu=reshape(velneuitrf(:,2)*1000,size(qx));
255+
qzv=qzv(ix,iy);
256+
qzu=qzu(ix,iy)./cosd(qyy);
257+
hold on
258+
quiverm(qyy(:),qxx(:),qzv(:),qzu(:),'r','filled')
259+
260+
%title(['Horizontal ' plttitle])
261+
title('Horizontal')
262+
263+
% Fig 4 - Plot the Vertical velocity/difference
264+
265+
qz=reshape(velneuitrf(:,3)*1000,size(qx));
266+
levels=getlevels(qz);
267+
268+
nexttile;
269+
%figure;
270+
ax=worldmap(latrange,lonrange);
271+
setm(ax,'GLineStyle',':','GColor','k','FFaceColor','none'); % 'FFaceColor','cyan'
272+
%[cs,h]=contourfm(qy,qx,qz,'LevelStep',levelstep,'LineColor',[0.6 0.6 0.6]);
273+
try
274+
[cs,h]=contourfm(qy,qx,qz,levels,'LineColor',[0.6 0.6 0.6]);
275+
ht = clabelm(cs,h);
276+
set(ht,'Color','white','BackgroundColor','none','FontWeight','normal')
277+
catch
278+
disp('Contourfm failure, skip...')
279+
end
280+
c=colorbar;
281+
c.Label.String=units;
282+
geoshow(landareas, 'FaceColor','none','EdgeColor','black');
283+
geoshow(rivers, 'Color', 'cyan')
284+
geoshow(borders, 'Color','black')
285+
%title(['Vertical ' plttitle])
286+
title('Vertical')
287+
288+
toc
289+
290+
end
291+
292+
function levels=getlevels(qz)
293+
294+
minz=min(qz(:));
295+
maxz=max(qz(:));
296+
zrange=maxz-minz;
297+
if zrange < 0.01
298+
minz=minz-0.01/2;
299+
maxz=maxz+0.01/2;
300+
zrange=maxz-minz;
301+
end
302+
levelstep=10^(floor(log10(zrange)));
303+
numlevels=ceil(zrange/levelstep);
304+
%levelstep=levelstep/floor(10/numlevels);
305+
while numlevels < 6
306+
levelstep=levelstep/2;
307+
numlevels=ceil(zrange/levelstep);
308+
end
309+
levelstep=max(levelstep,0.01);
310+
levels=floor(minz/levelstep)*levelstep:levelstep:ceil(maxz/levelstep)*levelstep;
311+
if length(levels) <1
312+
levels(1) = levels(1) - levelstep/2;
313+
levels(2) = levels(1) + levelstep;
314+
end
315+
316+
end

‎liveitrfmap.mlx‎

247 KB
Binary file not shown.

0 commit comments

Comments
 (0)