Patch with hole in map
<p>I would like to achieve the following result:</p>

<p><a href="https://i.stack.imgur.com/DmtdU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DmtdU.png" alt="Patch map"></a></p>

<p>As you can see, I would like to patch the entire map excluding the region formed by the union of the <a href="https://en.wikipedia.org/wiki/Circle_of_a_sphere" rel="nofollow noreferrer">small circles</a> that originate from location A and location B.</p>

<p>Locations A and B have the following coordinates (expressed in degrees):</p>

<pre><code>% Coordinates.
A = [43.6350000000000, 1.36777777777778];
B = [52.7019444444445, -8.92472222222222];
</code></pre>

<p>I obtain the coordinates of the small circles using the <a href="https://mathworks.com/help/map/ref/scircle1.html" rel="nofollow noreferrer"><code>scircle1</code></a> function. Both small circles have a radius of 654 nautical miles:</p>

<pre><code>% Small circles from center, range, and azimuth.
RangeNM = 654;
[latcA, loncA] = scircle1(A(1), A(2), RangeNM, [], earthRadius('nm'));
[latcB, loncB] = scircle1(B(1), B(2), RangeNM, [], earthRadius('nm'));
</code></pre>

<p>The <strong>internal contour</strong> of the patch is the perimeter formed by the coordinates of the <a href="https://en.wikipedia.org/wiki/Union_(set_theory)" rel="nofollow noreferrer">union</a> of both small circles. I obtain the coordinates using the <a href="https://mathworks.com/help/map/ref/polybool.html" rel="nofollow noreferrer"><code>polybool</code></a> function:</p>

<pre><code>% Set operations on polygonal regions.
[latU, lonU] = polybool('union', latcA, loncA, latcB, loncB);
</code></pre>

<p>The <strong>external contour</strong> of the patch is the perimeter formed by the latitude and longitude limits of the map. I use <a href="https://mathworks.com/help/map/ref/getm.html" rel="nofollow noreferrer"><code>getm</code></a> to obtain the <code>MapLatLimit</code> and <code>MapLonLimit</code> properties of the map. From the picture above, we can anticipate that <code>LatLim = [30, 70]</code> (30°N to 70°N) and <code>LonLim = [-30, 20]</code> (30°W to 20°E):</p>

<pre><code>% Get Map Latitude and Longitude limit.
LatLim = getm(gca, 'MapLatLimit');
LonLim = getm(gca, 'MapLonLimit');
</code></pre>

<p>Finally I try to create the patch using the <a href="https://mathworks.com/help/map/ref/patchm.html" rel="nofollow noreferrer"><code>patchm</code></a> function, which is the equivalent of the patch function for maps. Here is where I am having problems. I have tried three different approaches, but none of them is successfull:</p>

<pre><code>% APPROACH 1
lat = {[LatLim(1) LatLim(2) LatLim(2) LatLim(1) LatLim(1)], latU};
lon = {[LonLim(1) LonLim(1) LonLim(2) LonLim(2) LonLim(1)], lonU};
% Compute face and vertex matrices.
[f, v] = poly2fv(lat, lon);
patchm('Faces', f, 'Vertices', v, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4); %Doesn't work

% APPROACH 2
lat = {[LatLim(1) LatLim(2) LatLim(2) LatLim(1) LatLim(1)], latU};
lon = {[LonLim(1) LonLim(1) LonLim(2) LonLim(2) LonLim(1)], lonU};
% Compute face and vertex matrices.
[f, v] = poly2fv(lat, lon);
patchm(v(:,1), v(:,2), 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4); % Doesn't work

% APPROACH 3
lat = [latU', LatLim(1), LatLim(2), LatLim(2), LatLim(1), LatLim(1)];
lon = [lonU', LonLim(1), LonLim(1), LonLim(2), LonLim(2), LonLim(1)];
patchm(lat, lon, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4); % Doesn't work
</code></pre>

<p>I would really appreciate a helping hand.</p>

<hr>

<p>I have found a very good example of a patch with holes in the documentation page of the <a href="https://mathworks.com/help/map/ref/poly2fv.html" rel="nofollow noreferrer"><code>poly2fv</code></a> function. However, the example uses the standard <code>patch</code> function with cartesian coordinates (x, y). Note that <code>patch('Faces', f, 'Vertices', v,...</code> is used. I am trying to replicate this same example using the <code>patchm</code> function and geographic coordinates (lat, lon).</p>

<hr>

<h2>Code:</h2>

<blockquote>
  <p>Note: This code requires the <a href="https://mathworks.com/products/mapping/" rel="nofollow noreferrer">Mapping Toolbox</a>.</p>
</blockquote>

<pre><code>% Read vector features and attributes from shapefile.
landareas = shaperead('landareas.shp', 'UseGeoCoords', true);

% Define map axes and set map properties.
axesm ('lambert',...
    'MapLonLimit', [-30 20],...
    'MapLatLimit', [30 70],...
    'MapParallels', [38.00555556 71.01111111],...
    'Frame', 'on',...
    'FLineWidth', 1,...
    'Grid', 'on',...
    'GLineStyle', '-',...
    'GLineWidth', 0.1,...
    'GColor', [.7 .7 .7]);

% Display map latitude and longitude data.
geoshow(landareas, 'FaceColor', [1 1 .5], 'EdgeColor', [.3 .3 .3]);

% Toggle and control display of graticule lines.
gridm('MLineLocation', 5,...
    'MLabelLocation', 5,...
    'PLineLocation', 5,...
    'PLabelLocation', 5);

% Toggle and control display of meridian labels.
mlabel on;

% Toggle and control display of parallel labels.
plabel on;

axis off;

% Coordinates.
A = [43.6350000000000, 1.36777777777778];
B = [52.7019444444445, -8.92472222222222];

% Plot A.
plotm(A(1), A(2), '.r');
textm(A(1), A(2)+1, 'A');

% Plot B.
plotm(B(1), B(2), '.r');
textm(B(1), B(2)+1, 'B');

% Small circles from center, range, and azimuth.
RangeNM = 654;
[latcA, loncA] = scircle1(A(1), A(2), RangeNM, [], earthRadius('nm'));
[latcB, loncB] = scircle1(B(1), B(2), RangeNM, [], earthRadius('nm'));

% Set operations on polygonal regions.
[latU, lonU] = polybool('union', latcA, loncA, latcB, loncB);

% Get Map Latitude and Longitude limit.
LatLim = getm(gca, 'MapLatLimit');
LonLim = getm(gca, 'MapLonLimit');

% APPROACH 1
lat = {[LatLim(1) LatLim(2) LatLim(2) LatLim(1) LatLim(1)], latU};
lon = {[LonLim(1) LonLim(1) LonLim(2) LonLim(2) LonLim(1)], lonU};
% Compute face and vertex matrices.
[f, v] = poly2fv(lat, lon);
patchm('Faces', f, 'Vertices', v, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4); % Doesn't work

% APPROACH 2
lat = {[LatLim(1) LatLim(2) LatLim(2) LatLim(1) LatLim(1)], latU};
lon = {[LonLim(1) LonLim(1) LonLim(2) LonLim(2) LonLim(1)], lonU};
% Compute face and vertex matrices.
[f, v] = poly2fv(lat, lon);
patchm(v(:,1), v(:,2), 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4); % Doesn't work

% APPROACH 3
lat = [latU', LatLim(1), LatLim(2), LatLim(2), LatLim(1), LatLim(1)];
lon = [lonU', LonLim(1), LonLim(1), LonLim(2), LonLim(2), LonLim(1)];
patchm(lat, lon, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4); % Doesn't work
</code></pre>