mapprofile
Interpolate between waypoints on regular data grid
Syntax
[zi,ri,lat,lon] = mapprofile
[zi,ri,lat,lon] = mapprofile(Z,R,lat,lon)
[zi,ri,lat,lon] = mapprofile(Z,R,lat,lon,units)
[zi,ri,lat,lon] = mapprofile(Z,R,lat,lon,ellipsoid)
[zi,ri,lat,lon] = mapprofile(___,'trackmethod','interpmethod')
Description
mapprofile 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 grid's zdata is
used for the profile. The color data is used in the absence of zdata.
The result is displayed in a new figure.
[zi,ri,lat,lon] = mapprofile returns the
values of the profile without displaying them. The output zi contains
interpolated values along great circles between the waypoints. ri 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.
[zi,ri,lat,lon] = mapprofile(Z,R,lat,lon) accepts
as input a regular data grid and waypoint vectors. No displayed grid
is required. Sets of waypoints may be separated by NaNs into line
sequences. The output ranges are measured from the first waypoint
within a sequence. R can be a geographic raster
reference object, a referencing vector, or a referencing matrix.
If R is a geographic raster reference object,
its RasterSize property must be consistent with size(Z).
If R is a referencing vector, it must be
a 1-by-3 with elements:
[cells/degree northern_latitude_limit western_longitude_limit]
If R is a referencing matrix, it must be
3-by-2 and transform raster row and column indices to or from geographic
coordinates according to:
[lon lat] = [row col 1] * R
If R is a referencing matrix, it must define
a (non-rotational, non-skewed) relationship in which each column of
the data grid falls along a meridian and each row falls along a parallel.
[zi,ri,lat,lon] = mapprofile(Z,R,lat,lon,units) specifies
the units of the output ranges along the profile. Valid range units
inputs are any distance value recognized by unitsratio.
Surface distances are computed using the default radius of the earth.
If omitted, 'degrees' are assumed.
[zi,ri,lat,lon] = mapprofile(Z,R,lat,lon,ellipsoid)
uses the provided ellipsoid definition in computing the range along the profile.
ellipsoid is a referenceSphere, referenceEllipsoid, or oblateSpheroid object, or a vector of the form [semimajor_axis
eccentricity]. The output range is reported in the same distance units as the
semimajor axes of the ellipsoid. If omitted, the range vector is for a sphere.
[zi,ri,lat,lon] = mapprofile(___,'trackmethod','interpmethod') control
the interpolation methods used. Valid track methods are 'gc' for
great circle tracks between waypoints, and 'rh' for
rhumb lines. Valid methods for interpolation within the matrix 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 coastlines
plotm(coastlat,coastlon)
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 enables you to interactively
pick waypoints on the figure using the mouse. Each click selects a new waypoint and
stores the profile values, but the points are not displayed on the figure. Press
Enter after you select the final point:
axes(h); meshm(map, refvec, size(map)) demcmap(map) [zi,ri,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(ri),zi/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 = strcmp('Seattle', {cities(:).Name});
Frankfurt = strcmp('Frankfurt', {cities(:).Name});
lat = [cities(Seattle).Lat cities(Frankfurt).Lat]
lon = [cities(Seattle).Lon cities(Frankfurt).Lon]
load topo
[valp,rp,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