Documentation |
Intersection of two latitude-longitude quadrangles
[latlim, lonlim] = intersectgeoquad(latlim1, lonlim1,
latlim2, lonlim2)
[latlim, lonlim] = intersectgeoquad(latlim1, lonlim1, latlim2, lonlim2) computes the intersection of the quadrangle defined by the latitude and longitude limits latlim1 and lonlim1, with the quadrangle defined by the latitude and longitude limits latlim2 and lonlim2. latlim1 and latlim2 are two-element vectors of the form [southern-limit northern-limit]. Likewise, lonlim1 and lonlim2 are two-element vectors of the form [western-limit eastern-limit]. All input and output angles are in units of degrees. The intersection results are given in the output arrays latlim and lonlim. Given an arbitrary pair of input quadrangles, there are three possible results:
The quadrangles fail to intersect. In this case, both latlim and lonlim are empty arrays.
The intersection consists of a single quadrangle. In this case, latlim (like latlim1 and latlim2) is a two-element vector that has the form [southern-limit northern-limit], where southern-limit and northern-limit represent scalar values. lonlim (like lonlim1 and lonlim2), is a two-element vector that has the form [western-limit eastern-limit], with a pair of scalar limits.
The intersection consists of a pair of quadrangles. This can happen when longitudes wrap around such that the eastern end of one quadrangle overlaps the western end of the other and vice versa. For example, if lonlim1 = [-90 90] and lonlim2 = [45 -45], there are two intervals of overlap: [-90 -45] and [45 90]. These limits are returned in lonlim in separate rows, forming a 2-by-2 array. In our example (assuming that the latitude limits overlap), lonlim would equal [-90 -45; 45 90]. It still has the form [western-limit eastern-limit], but western-limit and eastern-limit are 2-by-1 rather than scalar. The two output quadrangles have the same latitude limits, but these are replicated so that latlim is also 2-by-2.
To continue the example, if latlim1 = [0 30] and latlim2 = [20 50], latlim equals [20 30; 20 30]. The form is still [southern-limit northern-limit], but in this case southern-limit and northern-limit are 2-by-1.
Nonintersecting quadrangles:
[latlim, lonlim] = intersectgeoquad( ... [-40 -60], [-180 180], [40 60], [-180 180]) latlim = [] lonlim = []
Intersection is a single quadrangle:
[latlim, lonlim] = intersectgeoquad( ... [-40 60], [-120 45], [-60 40], [160 -75]) latlim = -40 40 lonlim = -120 -75
Intersection is a pair of quadrangles:
[latlim, lonlim] = intersectgeoquad( ... [-30 90],[-10 -170],[-90 30],[170 10]) latlim = -30 30 -30 30 lonlim = -10 10 170 -170
Inputs and output fully encircle the planet:
[latlim, lonlim] = intersectgeoquad( ... [-30 90],[-180 180],[-90 30],[0 360]) latlim = -30 30 lonlim = -180 180
Find and map the intersection of the bounding boxes of adjoining U.S. states:
usamap({'Minnesota','Wisconsin'}) S = shaperead('usastatehi','UseGeoCoords',true,'Selector',... {@(name) any(strcmp(name,{'Minnesota','Wisconsin'})), 'Name'}); geoshow(S, 'FaceColor', 'y') textm([S.LabelLat], [S.LabelLon], {S.Name},... 'HorizontalAlignment', 'center') latlimMN = S(1).BoundingBox(:,2)' latlimMN = 43.4995 49.3844 lonlimMN = S(1).BoundingBox(:,1)' lonlimMN = -97.2385 -89.5612 latlimWI = S(2).BoundingBox(:,2)' latlimWI = 42.4918 47.0773 lonlimWI = S(2).BoundingBox(:,1)' lonlimWI = -92.8892 -86.8059 [latlim lonlim] = ... intersectgeoquad(latlimMN, lonlimMN, latlimWI, lonlimWI) latlim = 43.4995 47.0773 lonlim = -92.8892 -89.5612 geoshow(latlim([1 2 2 1 1]), lonlim([1 1 2 2 1]), ... 'DisplayType','polygon','FaceColor','m')