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')

`latlim1`

and `latlim2`

should
normally be given in order of increasing numerical value. No error
will result if, for example, `latlim1(2) < latlim1(1)`

,
but the outputs will both be empty arrays.

No such restriction applies to `lonlim1`

and `lonlim2`

.
The first element is always interpreted as the western limit, even
if it exceeds the second element (the eastern limit). Furthermore, `intersectgeoquad`

correctly
handles whatever longitude-wrapping convention may have been applied
to `lonlim1`

and `lonlim2`

.

In terms of output, `intersectgeoquad`

wraps `lonlim`

such
that all elements fall in the closed interval `[-180 180]`

.
This means that if (one of) the output quadrangle(s) crosses the 180°
meridian, its western limit exceeds its eastern limit. The result
would be such that

lonlim(2) < lonlim(1)

lonlim(k,2) < lonlim(k,1)

`k`

equals
1 or 2 if the intersection comprises a pair of quadrangles.If `abs(diff(lonlim1))`

or `abs(diff(lonlim2))`

equals `360`

,
its quadrangle is interpreted as a latitudinal zone that fully encircles
the planet, bounded only by one parallel on the south and another
parallel on the north. If two such quadrangles intersect, `lonlim`

is
set to `[-180 180]`

.

If you want to display geographic quadrangles generated by this
function or any other which are more than one or two degrees in extent,
they may not follow curved meridians and parallels very well. The
degree of departure depends on the extent of the quadrangle, the map
projection, and the map scale. In such cases, you can interpolate
intermediate vertices along quadrangle edges with the `outlinegeoquad`

function.