地图上有洞的补丁
Patch with hole in map
我想实现以下结果:
如您所见,我想修补整个地图,但不包括源自位置 A 和位置 B 的 small circles 的并集所形成的区域。
位置 A 和 B 的坐标如下(以度表示):
% Coordinates.
A = [43.6350000000000, 1.36777777777778];
B = [52.7019444444445, -8.92472222222222];
我使用scircle1
函数获取小圆圈的坐标。两个小圆的半径都是654海里:
% 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'));
patch的内部轮廓是由union of both small circles. I obtain the coordinates using the polybool
函数坐标形成的周长:
% Set operations on polygonal regions.
[latU, lonU] = polybool('union', latcA, loncA, latcB, loncB);
patch的外部轮廓是由地图的经纬度界限形成的周长。我使用 getm
获取地图的 MapLatLimit
和 MapLonLimit
属性。从上图中,我们可以预测 LatLim = [30, 70]
(30°N 至 70°N)和 LonLim = [-30, 20]
(30°W 至 20°E):
% Get Map Latitude and Longitude limit.
LatLim = getm(gca, 'MapLatLimit');
LonLim = getm(gca, 'MapLonLimit');
最后我尝试使用 patchm
函数创建补丁,这相当于地图的补丁函数。这是我遇到问题的地方。我尝试了三种不同的方法,但其中 none 是成功的:
% 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
非常感谢您的帮助。
我在 poly2fv
函数的文档页面中找到了一个很好的漏洞补丁示例。但是,该示例使用具有笛卡尔坐标 (x, y) 的标准 patch
函数。请注意,使用了 patch('Faces', f, 'Vertices', v,...
。我正在尝试使用 patchm
函数和地理坐标(纬度、经度)复制同一个示例。
代码:
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
似乎patchm
忽略了漏洞。一种解决方法是将周围的多边形分成两部分,以便区域的联合与两者重叠。然后从北部和南部补丁中减去该区域,绘制两个部分,并添加一个漂亮的边缘:
% 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);
我想实现以下结果:
如您所见,我想修补整个地图,但不包括源自位置 A 和位置 B 的 small circles 的并集所形成的区域。
位置 A 和 B 的坐标如下(以度表示):
% Coordinates.
A = [43.6350000000000, 1.36777777777778];
B = [52.7019444444445, -8.92472222222222];
我使用scircle1
函数获取小圆圈的坐标。两个小圆的半径都是654海里:
% 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'));
patch的内部轮廓是由union of both small circles. I obtain the coordinates using the polybool
函数坐标形成的周长:
% Set operations on polygonal regions.
[latU, lonU] = polybool('union', latcA, loncA, latcB, loncB);
patch的外部轮廓是由地图的经纬度界限形成的周长。我使用 getm
获取地图的 MapLatLimit
和 MapLonLimit
属性。从上图中,我们可以预测 LatLim = [30, 70]
(30°N 至 70°N)和 LonLim = [-30, 20]
(30°W 至 20°E):
% Get Map Latitude and Longitude limit.
LatLim = getm(gca, 'MapLatLimit');
LonLim = getm(gca, 'MapLonLimit');
最后我尝试使用 patchm
函数创建补丁,这相当于地图的补丁函数。这是我遇到问题的地方。我尝试了三种不同的方法,但其中 none 是成功的:
% 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
非常感谢您的帮助。
我在 poly2fv
函数的文档页面中找到了一个很好的漏洞补丁示例。但是,该示例使用具有笛卡尔坐标 (x, y) 的标准 patch
函数。请注意,使用了 patch('Faces', f, 'Vertices', v,...
。我正在尝试使用 patchm
函数和地理坐标(纬度、经度)复制同一个示例。
代码:
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
似乎patchm
忽略了漏洞。一种解决方法是将周围的多边形分成两部分,以便区域的联合与两者重叠。然后从北部和南部补丁中减去该区域,绘制两个部分,并添加一个漂亮的边缘:
% 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);