Remove Longitude Coordinate Discontinuities at Date Line Crossings

This example shows how to remove longitude coordinate discontinuities at date line crossings that can confuse set operations on polygons. This can happen when points with longitudes near 180 degrees connect to points with longitudes near -180 degrees, as might be the case for eastern Siberia and Antarctica, and also for small circles and other patch objects. To prepare geographic data for use with polybool or for patch rendering, cut the polygons at the date line with the flatearthpoly function. flatearthpoly returns a polygon with points inserted to follow the date line up to the pole, traverse the longitudes at the pole, and return to the date line crossing along the other edge of the date line.

Note: The toolbox display functions automatically cut and trim geographic data if required by the map projection. Use flatearthpoly only when performing set operations on polygons.

Create an orthographic view of the Earth and plot the coastlines on it.

axesm ortho
setm(gca,'Origin', [60 170]); framem on; gridm on
load coastlines
plotm(coastlat,coastlon)

Generate a small circle that encompasses the North Pole and color it yellow.

[latc,lonc] = scircle1(75,45,30);
patchm(latc,lonc,'y')

Flatten the small circle using the flatearthpoly function.

[latf,lonf] = flatearthpoly(latc,lonc);

Plot the cut circle that you just generated as a magenta line.

plotm(latf,lonf,'m')

Generate a second small circle that does not include a pole.

[latc1 lonc1] = scircle1(20, 170, 30);

Flatten the circle and plot it as a red line. Note that the second small circle, which does not cover a pole, is clipped into two pieces along the date line. The polygon for the first small circle is plotted in plane coordinates to illustrate its flattened shape. The flatearthpoly function assumes that the interior of the polygon being flattened is in the hemisphere that contains most of its edge points. Thus a polygon produced by flatearthpoly does not cover more than a hemisphere.

[latf1,lonf1] = flatearthpoly(latc1,lonc1);
plotm(latf1,lonf1,'r')