Search code examples
matlabplotmatlab-figuregeo

Patch with hole in map

I would like to achieve the following result:

Patch map

As you can see, I would like to patch the entire map excluding the region formed by the union of the small circles that originate from location A and location B.

Locations A and B have the following coordinates (expressed in degrees):

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

I obtain the coordinates of the small circles using the scircle1 function. Both small circles have a radius of 654 nautical miles:

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

The internal contour of the patch is the perimeter formed by the coordinates of the union of both small circles. I obtain the coordinates using the polybool function:

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

The external contour of the patch is the perimeter formed by the latitude and longitude limits of the map. I use getm to obtain the MapLatLimit and MapLonLimit properties of the map. From the picture above, we can anticipate that LatLim = [30, 70] (30°N to 70°N) and LonLim = [-30, 20] (30°W to 20°E):

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

Finally I try to create the patch using the patchm 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:

% 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

I would really appreciate a helping hand.

I have found a very good example of a patch with holes in the documentation page of the poly2fv function. However, the example uses the standard patch function with cartesian coordinates (x, y). Note that patch('Faces', f, 'Vertices', v,... is used. I am trying to replicate this same example using the patchm function and geographic coordinates (lat, lon).

Code:

Note: This code requires the Mapping Toolbox.

% 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
Solution

  • It seems like patchm ignores holes. A workaround would be to split the surrounding polygon into two parts, so that the union of regions overlaps with both. Then subtract that region from the northern and the southern patch, plot the two parts, and add a nice edge:

    % Find a point where we can split the polygon
split_lat = median(latU);
% generate northern patch
p_up_lt = [split_lat LatLim([2 2]) split_lat];
p_up_lg = LonLim([1 1 2 2]);

[lonAup, latAup] = polybool('-', p_up_lg, p_up_lt, lonU, latU);
% and plot it, with invisible edge
p_up = patchm(latAup, lonAup, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4, 'EdgeColor', 'none');

% generate and plot southern patch
p_dn_lt =  [LatLim(1) split_lat split_lat LatLim(1)];
p_dn_lg =  LonLim([1 1 2 2]);
[lonAdn, latAdn] = polybool('-', p_dn_lg, p_dn_lt, lonU, latU);
p_dn = patchm(latAdn, lonAdn, 'FaceColor', [.5 .5 .5], 'FaceAlpha', 0.4, 'EdgeColor', 'none');

% plot edge for the union
p_e = patchm(latU, lonU, 'FaceColor', 'none', 'FaceAlpha', 0.4);