Mapping Toolbox™

mapprofile - Interpolate heights between waypoints on regular data grid

Syntax

[z,rng,lat,lon] = mapprofile [z,rng,lat,lon] = mapprofile(Z,refvec,lat,lon) [z,rng,lat,lon] = mapprofile(Z,refvec,lat,lon,rngunits) [z,rng,lat,lon] = mapprofile(Z,refvec,lat,lon,ellipsoid) [z,rng,lat,lon] = ... mapprofile(Z,refvec,lat,lon,rngunits,trackmethod,interpmethod) [z,rng,lat,lon] = ... mapprofile(Z,refvec,lat,lon,ellipsoid,trackmethod,interpmethod)

Description

mapprofile computes or plots a profile of values between waypoints on a displayed regular data grid. mapprofile uses the current object if it is a regular data grid, or the first regular data grid found on the current axes. The map's ZData is used for the profile. The color data is used in the absence of data in z. The result is displayed in a new 3-D figure. When no input arguments are given, cross hairs appear on the map in the current figure. Digitize a transect with as many waypoints as you need, and after the final one press Enter to complete the command.

[z,rng,lat,lon] = mapprofile returns the values of the profile without displaying them. The output z contains interpolated values from ZData along great circles between the waypoints. rng is a vector of associated distances from the first waypoint in units of degrees of arc along the surface. lat and lon are the corresponding latitudes and longitudes. This too is an interactive version of mapprofile, in which you digitize a transect on the current map.

[z,rng,lat,lon] = mapprofile(Z,refvec,lat,lon) takes as input a regular data grid, Z, a three-element referencing vector, refvec, and waypoint vectors. No displayed map is required. Sets of waypoints can be separated by NaNs into line sequences. The output ranges are measured from the first waypoint within a sequence.

[z,rng,lat,lon] = mapprofile(Z,refvec,lat,lon,rngunits) specifies the units of the output ranges along the profile. Valid range units inputs are any distance string recognized by unitsratio. Surface distances are computed using the default radius of the grid. If omitted, 'degrees' is assumed.

[z,rng,lat,lon] = mapprofile(Z,refvec,lat,lon,ellipsoid) uses the provided ellipsoid definition in computing the range along the profile. The ellipsoid vector is of the form [semimajor axes, eccentricity]. The output range is reported in the same distance units as the semimajor axes of the ellipsoid vector. If omitted, the range vector is for a sphere.

[z,rng,lat,lon] = ... mapprofile(Z,refvec,lat,lon,rngunits,trackmethod,interpmethod) and [z,rng,lat,lon] = ... mapprofile(Z,refvec,lat,lon,ellipsoid,trackmethod,interpmethod) control the interpolation methods used. Valid trackmethods are 'gc' for great circle tracks between waypoints, and 'rh' for rhumb lines. Valid interpmethods for interpolation within the data grid are 'bilinear' for linear interpolation, 'bicubic' for cubic interpolation, and 'nearest' for nearest neighbor interpolation. If omitted, 'gc' and 'bilinear' are assumed.

Examples

Example 1

Create a map axes for the Korean peninsula. Specify an elevation profile across the sample Korean digital elevation data and plot it, combined with a coastline and city markers:

load korea
h = worldmap(map, refvec); % The figure has no map content
plat = [ 43  43  41  38];
plon = [116 120 126 128];
mapprofile(map, refvec, plat, plon)
load coast
plotm(lat, long)
geoshow('worldcities.shp', 'Marker', '.', 'Color', 'red')

When you select more than two waypoints, the automatically generated figure displays the result in three dimensions. The following example shows the relative sizes of the mountains in northern China compared to the depths of the Sea of Japan. The call to mapprofile without input arguments requires you to interactively pick waypoints on the figure using the mouse, and press Enter after you select the final point:

axes(h);
meshm(map, refvec, size(map))
demcmap(map)
[z,rng,lat,lon] = mapprofile;

Adding output arguments suppresses the display of the results in a new figure. You can then use the results in further calculations or display the results yourself. Here the profile from the upper left to lower right is computed from waypoints interactively picked on the map (your profile will not be identical to what is shown below). The example converts ranges and elevations to kilometers and displays them in a new figure, setting the vertical exaggeration factor to 20. With no vertical exaggeration, the changes in elevation would be almost too small to see.

figure
plot(deg2km(rng),z/1000)
daspect([ 1 1/20 1 ]);
xlabel 'Range (km)'
ylabel 'Elevation (km)'

Naturally, the profile you get depends on the transect locations you pick.

Example 2

You can compute values along a path without reference to an existing figure by providing a regular data grid and vectors of waypoint coordinates. Optional arguments allow control over the units of the range output and interpolation methods between waypoints and data grid elements.

Show what land and ocean areas lie under a great circle track from Frankfurt to Seattle:

cities = shaperead('worldcities.shp', 'UseGeoCoords', true);
Seattle = strmatch('Seattle', {cities(:).Name});
Frankfurt = strmatch('Frankfurt', {cities(:).Name});
lat = [cities(Seattle).Lat cities(Frankfurt).Lat]
lon = [cities(Seattle).Lon cities(Frankfurt).Lon]
load topo
[valp,rngp,latp,lonp] = ...
   mapprofile(double(topo),topolegend, ...
              lat,lon,'km','gc','nearest');
figure
worldmap([40 80],[-135 20])
land = shaperead('landareas.shp', 'UseGeoCoords', true);
faceColors = makesymbolspec('Polygon',...
   {'INDEX', [1 numel(land)], 'FaceColor', ...
   polcmap(numel(land))});
geoshow(land,'SymbolSpec',faceColors)
plotm(latp,lonp,'r')
plotm(lat,lon,'ro')
axis off