Documentation Center |
On this page… |
---|
Geographic Interpolation of Vectors Overlaying Polygons with Set Logic Cutting Polygons at the Date Line Trimming Vector Data to a Rectangular Region |
It can be difficult to identify line or patch segments once they have been combined into large NaN-clipped vectors. You can separate these polygon or line vectors into their component segments using polysplit, which takes column vectors as inputs.
Enter two NaN-delimited arrays in the form of column vectors:
lat = [45.6 -23.47 78 NaN 43.9 -67.14 90 -89]'; long = [13 -97.45 165 NaN 0 -114.2 -18 0]';
Use polysplit to create two cell arrays, latc and lonc:
[latc,lonc] = polysplit(lat,long) latc = [3x1 double] [4x1 double] lonc = [3x1 double] [4x1 double]
Inspect the contents of the cell arrays:
[latc{1} lonc{1}] [latc{2} lonc{2}] ans = 45.6 13 -23.47 -97.45 78 165 ans = 43.9 0 -67.14 -114.2 90 -18 -89 0
Note that each cell array element contains a segment of the original line.
To reverse the process, use polyjoin:
[lat2,lon2] = polyjoin(latc,lonc);
The joined segments are identical with the initial lat and lon arrays:
[lat long] == [lat2 lon2] ans = 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1
The logical comparison is false for the NaN delimiters, by definition.
You can test for global equality, including NaNs, as follows:
isequaln(lat,lat2) & isequaln(long,lon2) ans = 1
See the polysplit and polyjoin reference pages for further information.
A common operation on sets of line segments is the concatenation of segments that have matching endpoints. The polymerge command compares endpoints of segments within latitude and longitude vectors to identify endpoints that match exactly or lie within a specified distance. The matching segments are then concatenated, and the process continues until no more coincidental endpoints can be found. The two required arguments are a latitude (x) vector and a longitude (y) vector. The following exercise shows this process at work.
Construct column vectors representing coordinate values:
lat = [3 2 NaN 1 2 NaN 5 6 NaN 3 4]'; lon = [13 12 NaN 11 12 NaN 15 16 NaN 13 14]';
Concatenate the segments that match exactly:
[latm,lonm] = polymerge(lat,lon); [latm lonm] ans = 1 11 2 12 3 13 4 14 NaN NaN 5 15 6 16 NaN NaN
The original four segments are merged into two segments.
The polymerge function takes an optional third argument, a (circular) distance tolerance that permits inexact matching. A fourth argument enables you to specify whether the function outputs vectors or cell arrays. See the polymerge reference page for further information.
When using vector data, remember that, like raster data, coordinates are sampled measurements. This involves unavoidable assumptions concerning what the geographic reality is between specified data points. The normal assumption when plotting vector data requires that points be connected with straight line segments, which essentially indicates a lack of knowledge about conditions between the measured points. For lines that are by nature continuous, such as most rivers and coastlines, such piecewise linear interpolation can be false and misleading, as the following figure depicts.
Interpolating Sparse Vector Data
Despite the possibility of misinterpretation, circumstances do exist in which geographic data interpolation is useful or even necessary. To do this, use the interpm function to interpolate between known data points. One value of linearly interpolating points is to fill in lines of constant latitude or longitude (e.g., administrative boundaries) that can curve when projected.
This example interpolates values in a set of latitude and longitude points to have no more than one degree of separation in either direction.
Define two fictitious latitude and longitude data vectors:
lats = [1 2 4 5]; longs = [1 3 4 5];
Specify a densification parameter of 1 (the default unit is degrees):
maxdiff = 1;
Call interpm to fill in any gaps greater than 1º in either direction:
[newlats,newlongs] = interpm(lats,longs,maxdiff) newlats = 1.0000 1.5000 2.0000 3.0000 4.0000 5.0000 newlongs = 1.0000 2.0000 3.0000 3.5000 4.0000 5.0000
In lats, a gap of 2º exists between the values 2 and 4. A linearly interpolated point, (3,3.5) was therefore inserted in newlats and newlongs. Similarly, in longs, a gap of 2º exists between the 1 and the 3. The point (1.5, 2) was therefore interpolated and placed into newlats and newlongs. Now, the separation of adjacent points is no greater than maxdiff in either newlats or newlongs.
See the interpm reference page for further information.
interpm returns both the original data and new linearly interpolated points. Sometimes, however, you might want only the interpolated values. The functions intrplat and intrplon work similarly to the MATLAB^{®} interp1 function, and give you control over the method used for interpolation. Note that they only interpolate and return one value at a time.
Use intrplat to interpolate a latitude for a given longitude. Given a monotonic set of longitudes and their matching latitude points, you can interpolate a new latitude for a longitude you specify, interpolating along linear, spline, cubic, rhumb line, or great circle paths. The longitudes must increase or decrease monotonically. If this is not the case, you might be able to use the intrplon companion function if the latitude values are monotonic.
Interpolate a latitude corresponding to a longitude of 7.3º in the following data in a linear, great circle, and rhumb line sense:
Define some fictitious latitudes and longitudes:
longs = [1 3 4 9 13]; lats = [57 68 60 65 56];
Specify the longitude for which to compute a latitude:
newlong = 7.3;
Generate a new latitude with linear interpolation:
newlat = intrplat(longs,lats,newlong,'linear')
newlat =
63.3000
Generate the latitude using great circle interpolation:
newlat = intrplat(longs,lats,newlong,'gc')
newlat =
63.5029
Generate it again, specifying interpolation along a rhumb line:
newlat = intrplat(longs,lats,newlong,'rh')
newlat =
63.3937
The following diagram illustrates these three types of interpolation. The intrplat function also can perform spline and cubic spline interpolations.
As mentioned above, the intrplon function provides the capability to interpolate new longitudes from a given set of longitudes and monotonic latitudes.
See the intrplat and intrplon reference pages for further information.
A set of Mapping Toolbox™ functions perform intersection calculations on vector data computed by the toolbox, which include great and small circles as well as rhumb line tracks. The functions also determine intersections of arbitrary vector data.
Compute the intersection of a small circle centered at (0º,0º) with a radius of 1250 nautical miles and a small circle centered at (5ºN,30ºE) with a radius of 2500 kilometers:
[lat,long] = scxsc(0,0,nm2deg(1250),5,30,km2deg(2500)) lat = 17.7487 -12.9839 long = 11.0624 16.4170
In general, small circles intersect twice or never. For the case of exact tangency, scxsc returns two identical intersection points. Other similar commands include rhxrh for intersecting rhumb lines, gcxgc for intersecting great circles, and gcxsc for intersecting a great circle with a small circle.
Imagine a ship setting sail from Norfolk, Virginia (37ºN,76ºW), maintaining a steady due-east course (90º), and another ship setting sail from Dakar, Senegal (15ºN,17ºW), with a steady northwest course (315º). Where would the tracks of the two vessels cross?
[lat,long] = rhxrh(37,-76,90,15,-17,315) lat = 37 long = -41.7028
The intersection of the tracks is at (37ºN,41.7ºW), which is roughly 600 nautical miles west of the Azores in the Atlantic Ocean.
You can also compute the intersection points of arbitrary vectors of latitude and longitude. The polyxpoly command finds the segments that intersect and interpolates to find the intersection points. The interpolation is done linearly, as if the points were in a Cartesian x-y coordinate system. The polyxpoly command can also identify the line segment numbers associated with the intersections:
[xint,yint] = polyxpoly(x1,y1,x2,y2);
If the spacing between points is large, there can be some difference between the intersection points computed by polyxpoly and the intersections shown on a map display. This is a result of the difference between straight lines in the unprojected and projected coordinates. Similarly, there can be differences between the polyxpoly result and intersections assuming great circles or rhumb lines between points.
Use the function areaint to calculate geographic areas for vector data in polygon format. The function performs a numerical integration using Green's Theorem for the area on a surface enclosed by a polygon. Because this is a discrete integration on discrete data, the results are not exact. Nevertheless, the method provides the best means of calculating the areas of arbitrarily shaped regions. Better measures result from better data.
The Mapping Toolbox function areaint (for area by integration), like the other area functions, areaquad and areamat, returns areas as a fraction of the entire planet's surface, unless a radius is provided. Here you calculate the area of the continental United States using the conus MAT-file. Three areas are returned, because the data contains three polygons: Long Island, Martha's Vineyard, and the rest of the continental U.S.:
load conus earthradius = earthRadius('km'); area = areaint(uslat,uslon,earthradius) area = 1.0e+06 * 7.9256 0.0035 0.0004
Because the default Earth radius is in kilometers, the area is in square kilometers. From the same variables, the areas of the Great Lakes can be calculated, this time in square miles:
load conus earthradius = earthRadius('miles'); area = areaint(gtlakelat,gtlakelon,earthradius) area = 1.0e+004 * 8.0119 1.0381 0.7634
Again, three areas are returned, the largest for the polygon representing Superior, Michigan, and Huron together, the other two for Erie and Ontario.
Polygon set operations are used to answer a variety of questions about logical relationships of vector data polygon objects. Standard set operations include intersection, union, subtraction, and an exclusive OR operation. The polybool function performs these operations on two sets of vectors, which can represent x-y or latitude-longitude coordinate pairs. In computing points where boundaries intersect, interpolations are carried out on the coordinates as if they were planar. Here is an example that shows all the available operations.
The result is returned as NaN-clipped vectors by default. In cases where it is important to distinguish outer contours of polygons from interior holes, polybool can also accept inputs and return outputs as cell arrays. In the cell array format, a cell array entry starts with the list of points making up the outer contour. Subsequent NaN-clipped faces within the cell entry are interpreted as interior holes.
The following exercise demonstrates how you can use polybool:
Construct a twelve-sided polygon:
theta = -(0:pi/6:2*pi)'; lat1 = sin(theta); lon1 = cos(theta);
Construct a triangle that overlaps it:
lat2 = [0 1 -1 0]'; lon2 = [0 2 2 0]';
Plot the two shapes together with blue and red lines:
axesm miller plotm(lat1,lon1,'b') plotm(lat2,lon2,'r')
Compute the intersection polygon and plot it as a green patch:
[loni,lati] = polybool('intersection',lon1,lat1,lon2,lat2); [lati loni] geoshow(lati,loni,'DisplayType','polygon','FaceColor','g') ans = 0 1.0000 -0.4409 0.8819 0 0 0.4409 0.8819 0 1.0000
Compute the union polygon and plot it as a magenta patch:
[lonu,latu] = polybool('union',lon1,lat1,lon2,lat2); [latu lonu] geoshow(latu,lonu,'DisplayType','polygon','FaceColor','m') ans = -1.0000 2.0000 -0.4409 0.8819 -0.5000 0.8660 -0.8660 0.5000 -1.0000 0.0000 -0.8660 -0.5000 -0.5000 -0.8660 0 -1.0000 0.5000 -0.8660 0.8660 -0.5000 1.0000 -0.0000 0.8660 0.5000 0.5000 0.8660 0.4409 0.8819 1.0000 2.0000 -1.0000 2.0000
Compute the exclusive OR polygon and plot it as a yellow patch:
[lonx,latx] = polybool('xor',lon1,lat1,lon2,lat2); [latx lonx] geoshow(latx,lonx,'DisplayType','polygon','FaceColor','y') ans = -1.0000 2.0000 -0.4409 0.8819 -0.5000 0.8660 -0.8660 0.5000 -1.0000 0.0000 -0.8660 -0.5000 -0.5000 -0.8660 0 -1.0000 0.5000 -0.8660 0.8660 -0.5000 1.0000 -0.0000 0.8660 0.5000 0.5000 0.8660 0.4409 0.8819 1.0000 2.0000 -1.0000 2.0000 NaN NaN 0.4409 0.8819 0 0 -0.4409 0.8819 0 1.0000 0.4409 0.8819
Subtract the triangle from the 12-sided polygon and plot the resulting concave polygon as a white patch:
[lonm,latm] = polybool('minus',lon1,lat1,lon2,lat2); [latm lonm] geoshow(latm,lonm,'DisplayType','polygon','FaceColor','w') ans = 0.8660 0.5000 0.5000 0.8660 0.4409 0.8819 0 0 -0.4409 0.8819 -0.5000 0.8660 -0.8660 0.5000 -1.0000 0.0000 -0.8660 -0.5000 -0.5000 -0.8660 0 -1.0000 0.5000 -0.8660 0.8660 -0.5000 1.0000 -0.0000 0.8660 0.5000
The final set of colored shapes is shown below.
See the polybool reference page for further information.
Polygon set operations treat input vectors as plane coordinates. The polyxpoly function can be confused by geographic data that has discontinuities in longitude coordinates at date-line crossings. This can happen when points with longitudes near 180º connect to points with longitudes near -180º, as might be the case for eastern Siberia and Antarctica, and also for small circles and other patch objects generated by toolbox functions.
You can prepare such geographic data for use with polybool or for patch rendering by cutting the polygons at the date line with the flatearthpoly function. The result of flatearthpoly is a polygon with points inserted to follow the date line up to the pole, traverse the longitudes at the pole, and return to the date line crossing along the other edge of the date line.
Create an orthographic view of the Earth and plot coast on it:
axesm ortho setm(gca,'Origin', [60 170]); framem on; gridm on load coast plotm(lat, long)
Generate a small circle that encompasses the North Pole and color it yellow:
[latc,lonc] = scircle1(75,45,30);
patchm(latc,lonc,'y')
Flatten the small circle with flatearthpoly:
[latf,lonf] = flatearthpoly(latc,lonc);
Plot the cut circle that you just generated as a magenta line:
plotm(latf,lonf,'m')
Generate a second small circle that does not include a pole:
[latc1 lonc1] = scircle1(20, 170, 30);
Flatten it and plot it as a red line:
[latf1,lonf1] = flatearthpoly(latc1,lonc1);
plotm(latf1,lonf1,'r')
The following figure shows the result of these operations. Note that the second small circle, which does not cover a pole, has been clipped into two pieces along the date line. On the right, the polygon for the first small circle is plotted in plane coordinates to illustrate its flattened shape.
The flatearthpoly function assumes that the interior of the polygon being flattened is in the hemisphere that contains most of its edge points. Thus a polygon produced by flatearthpoly does not cover more than a hemisphere.
Note As this figure illustrates, you do not need to use flatearthpoly to prepare data for a map display. The Mapping Toolbox display functions automatically cut and trim geographic data if required by the map projection. Use this function only when conducting set operations on polygons. |
See the flatearthpoly reference page for further information.
A buffer zone is the area within a specified distance of a map feature. For vector geodata, buffer zones are constructed as polygons. For raster geodata, buffer zones are collections of contiguous, identically coded grid cells. When the feature is a polygon, a buffer zone can be defined as the locus of points within a certain distance of its boundary, either inside or outside the polygon. A buffer zone for a linear object is the locus of points a certain distance away from it. Buffer zones form equidistant contour lines around objects.
The bufferm function computes and returns vectors that represent a set of points that define a buffer zone. It forms the buffer by placing small circles at the vertices of the polygon and rectangles along each of its line segments, and applying a polygon union set operation to these objects.
Demonstrate bufferm using a polygon representing the Island of Madagascar that you extract from the landareas data set. The boundary of Madagascar is passed to bufferm as latitude and longitude vectors. Using this data, compute a buffer zone at a distance of 0.75 degrees in from the boundaries of Madagascar:
Create a base map of the area surrounding Madagascar:
ax = worldmap('madagascar'); madagascar = shaperead('landareas',... 'UseGeoCoords', true,... 'Selector', {@(name)strcmpi(name,'Madagascar'), 'Name'}); geoshow(ax, madagascar)
Use bufferm to process the polygon and output a buffer zone .75 degrees inland:
[latb,lonb] = bufferm(madagascar.Lat, madagascar.Lon, .75, 'in');
Show the buffer zone in green:
geoshow(latb, lonb, 'DisplayType', 'polygon', 'FaceColor', 'green')
It is not unusual for vector data to extend beyond the geographic region currently of interest. For example, you might have coastline data for the entire world (such as the coast data set), but are interested in mapping Australia only. In this and other situations, you might want to eliminate unnecessary data from the workspace and from calculations in order to save memory or to speed up processing and display.
Line data and patch data need to be trimmed differently. You can trim line data by simply removing points outside the region of interest by clipping lines at the map frame or to some other defined region. Patch data requires a more complicated method to ensure that the patch objects are correctly formed.
For the vector data, two functions are available to achieve this. If the vectors are to be handled as line data, the maptriml function returns variables containing only those points that lie within the defined region. If, instead, you want to maintain polygon format, use the maptrimp function. Be aware, however, that patch-trimmed data is usually larger and more expensive to compute.
Note When drawing maps, Mapping Toolbox display functions automatically trim vector geodata to the region specified by the frame limits (FLatLimit and FLonLimit map axes properties) for azimuthal projections, or to frame or map limits (MapLatLimit and MapLonLimit map axes properties) for nonazimuthal projections. The trimming is done internally in the display routine, keeping the original data intact. For further information on trimming vector geodata, see Axes for Drawing Maps, along with the reference pages for the trimming functions. |
Load the coast MAT-file for the entire world:
load coast
Define a region of interest centered on Australia:
latlim = [-50 0]; longlim = [105 160];
Use maptriml to delete all data outside these limits, producing line vectors:
[linelat,linelong] = maptriml(lat,long,latlim,longlim);
Do this again, but use maptrimp to produce polygon vectors:
[polylat,polylong] = maptrimp(lat,long,latlim,longlim);
See how much data has been reduced:
whos
Name Size Bytes Class lat 9589x1 76712 double latlim 1x2 16 double linelat 870x1 6960 double linelong 870x1 6960 double long 9589x1 76712 double longlim 1x2 16 double polylat 878x1 7024 double polylong 878x1 7024 double
Note that the clipped data is only 10% as large as the original data set.
Plot the trimmed patch vectors on a Miller projection:
axesm('MapProjection', 'miller', 'Frame', 'on',... 'FlatLimit', latlim, 'FlonLimit', longlim) patchesm(polylat, polylong, 'c')
Plot the trimmed line vectors to see that they conform to the patches:
plotm(linelat, linelong, 'm')
See the maptriml and maptrimp reference pages for further information.
Often a set of data contains unwanted data mixed in with the desired values. For example, your data might include vectors covering the entire United States, but you only want to work with those falling in Alabama. Sometimes a data set contains noise—perhaps three or four points out of several thousand are obvious errors (for example, one of your city points is in the middle of the ocean). In such cases, locating outliers and errors in the data arrays can be quite tedious.
The filterm command uses a data grid to filter a vector data set. Its calling sequence is as follows:
[flats,flons] = filterm(lats,lons,grid,refvector,allowed)
Each location defined by lats and lons is mapped to a cell in grid, and the value of that grid cell is obtained. If that value is found in allowed, that point is output to flats and flons. Otherwise, the point is filtered out.
The grid might encode political units, and the allowed values might be the code or codes indexing certain states or countries (e.g., Alabama). The grid might also be real-valued (e.g., terrain elevations), although it could be awkward to specify all the values allowed. More often, logical or relational operators give better results for such grids, enabling the allowed value to be 1 (for true). For example, you could use this transformation of the topo grid:
[flats,flons] = filterm(lats,lons,double(topo>0),topolegend,1)
The output would be those points in lats and lons that occupy dry land (mostly because some water bodies are above sea level).
For further information, see the filterm reference page. Also see Data Grids as Logical Variables.
Avoiding visual clutter in composing maps is an essential part of cartographic presentation. In cartography, this is described as map generalization, which involves coordinating many techniques, both manual and automated. Limiting the number of points in vector geodata is an important part of generalizing maps, and is especially useful for conditioning cartographic data, plotting maps at small scales, and creating versions of geodata for use at small scales.
An easy, but naive, approach to point reduction is to discard every nth element in each coordinate vector (simple decimation). However, this can result in poor representations of the original shapes. The toolbox provides a function to eliminate insignificant geometric detail in linear and polygonal objects, while still maintaining accurate representations of their shapes. The reducem function implements a powerful line simplification algorithm (known as Douglas-Peucker) that intelligently selects and deletes visually redundant points.
The reducem function takes latitude and longitude vectors, plus an optional linear tolerance parameter as arguments, and outputs reduced (simplified) versions of the vectors, in which deviations perpendicular to local "trend lines" in the vectors are all greater than the tolerance criterion. Endpoints of vectors are preserved. Optional outputs are an error measure and the tolerance value used (it is computed when you do not supply a value).
Note Simplified line data might not always be appropriate for display. If all or most intermediate points in a feature are deleted, then lines that appear straight in one projection can be incorrectly displayed as straight lines in others, and separate lines can be caused to intersect. In addition, when you are reducing data over large world regions, the effective degree of reduction near the poles are less than that achieved near the equator, due to the fact that the algorithm treats geographic coordinates as if they were planar. |
The reducem function works on both patch and line data. Getting results that look right at an intended scale might require some experimentation with the tolerance parameter. The best way to proceed might be to allow the tolerance to default, and have reducem return the tolerance it computed as the fourth return argument. If the output still has too much detail, then double the tolerance and try again. Similarly, if the output lines do not have enough detail, halve the tolerance and try again. You can also use the third return parameter, which indicates the percentage of line length that was eliminated by reduction, as a guide to achieve consistent simplification results, although this parameter is sensitive to line geometry and thus can vary by feature type.
To demonstrate the use of reducem, this example extracts the outline of the state of Massachusetts from the usastatehi high-resolution shapefile:
Read Massachusetts data from the shapefile. Use the Selector parameter to read only the vectors representing the Massachusetts state boundaries:
ma = shaperead('usastatehi.shp',... 'UseGeoCoords', true,... 'Selector', {@(name)strcmpi(name,'Massachusetts'), 'Name'});
Extract the coordinate data for simplification. There are 957 points to begin with:
maLat = ma.Lat; maLon = ma.Lon; numel(maLat)
ans = 957
Use reducem to simplify the boundary vectors, and inspect the results:
[maLat1, maLon1, cerr, tol] = reducem(maLat', maLon'); numel(maLat1)
ans = 252
The number of points used to represent the boundary has dropped from 958 to 253. Compute the degree of reduction:
numel(maLat1)/numel(maLat)
ans = 0.2633
The vectors have been reduced to about a quarter of their original size using the default tolerance.
Examine the error and tolerance values returned by reducem:
[cerr tol]
ans = 0.0331 0.0060
The cerr value says that only 3.3% of total boundary length was eliminated (despite removing 74% of the points). The tolerance that achieved this was computed by reducem as 0.006 decimal degrees, or about 0.66 km.
Use geoshow to plot the reduced outline in red over the original outline in blue:
figure axesm('MapProjection', 'eqdcyl', 'MapLatLim', [41.1 43.0],... 'MapLonLim', [-73.6, -69.8], 'Frame', 'off', 'Grid', 'off'); geoshow(maLat, maLon, 'DisplayType', 'line', 'color', 'blue') geoshow(maLat1, maLon1, 'DisplayType', 'line', 'color', 'red')
Differences in details are not apparent unless you zoom in two or three times; click the Zoom tool to explore the map.
Double the tolerance, and reduce the original boundary into new variables:
[maLat2,maLon2,cerr2,tol2] = reducem(maLat', maLon', 0.012);
Repeat step 3 with new data and plot it in dark green:
numel(maLat2)/numel(maLat)
ans = 0.1641
[cerr2 tol2]
ans = 0.0517 0.0120
geoshow(maLat2, maLon2, 'DisplayType', 'line', 'color', [0 .6 0])
Now you have removed 84% of the points, and 5.2% of total length.
Repeat one more time, raising the tolerance to 0.1 degrees, and plot the result in black:
[maLat3, maLon3, cerr3, tol3] = reducem(maLat', maLon', 0.1); geoshow(maLat3, maLon3, 'DisplayType', 'line', 'color', 'black')
As overlaid with the original data, the reduced state boundaries look like this.
As this example and the composite map below demonstrate, the visual effects of point reduction are subtle, up to a point. Most of the vertices can be eliminated before the effects of line simplification are very noticeable. Experimentation is usually required, because the visual effects a given value for a tolerance yield depend on the degrees and types of line complexity, and they are often nonlinear with respect to tolerance values. Also, the extent of line detail reduction should be informed by the purpose of the map and the scale at which it is to be displayed.
Note This exercise generalized a set of disconnected patches. When patches are contiguous (such as the U.S. state outlines), using reducem can result in inconsistencies in boundary representation and gaps at points where states meet. For best results, reducem should be applied to line data. |
See the reducem reference page for further information.