polyshape fractures after multiple merges

I am collecting satellite weather data and plotting it on a Mercator projection Earth map. I need to determine how long it takes to cover the entire sea area of the Earth with these scan patches. I create a polyshape for each satellite scan, and then merge all of the polys into a single, overall polyshape which I compare against the "coastlines" poly, comparing the number of holes in each. When the hole numbers match, I assume that I have covered the whole Earth (land or water). Here is a depiction of my scan with a single scan polyshape:

When the scan is complete, I will have filled in the entire water area of the Earth. Any holes in the polyshape will represent continents and islands. The "coastlines" polyshape has 221 regions and 22 holes. My merged polyshape should match this.

Here is a simplified version of my code:

pgon_merge = polyshape;

gx = axes;
load coastlines;
coast = polyshape(coastlon,coastlat);
plot(gx,coast);

srcfiles = dir(filepath);

for i = 1,numl(srcFiles)
    
    % read lat, lon from source files directory, 1 file at a time
    % lat, lon size: (286x636)
    
    filename = srcFiles(i).name;
    
    lat = ncread(filename,'Latitude');
    lon = ncread(filename,'Longitude');
    
    k = boundary(lat(:),lon(:));
    pgon = polyshape(lon(k),lat(k));
    
    pgon_land = intersect(coast,pgon);	% land only patch
    pgon_sea = xor(pgon,pgon_land);		% sea only patch
    
    pgon_merge = union(pgon_sea,pgon);	% merge sea poly into global poly
    plot(gx,pgon_merge);
    
    disp(strcat("patch holes = ",num2str(pgon_merge.NumHoles)));
    
end

In order to plot just the land areas, I do an intersect on the coast poly and the latest scan poly. Conversely, to keep just the water areas, I do a "xor" on the latest scan poly and the land areas. Then I do a "union" to combine that poly with the growing merge_poly.

This works, up to a point, but after several hundred merges I start to get strange lines sweeping across the merged polyshape. I suspect these are artifacts of the merge process, where lat/lon points from one poly are combined with points from the merged poly in unexpected ways.

(Note: this is a hand-drawn depiction of what I'm seeing, not the actual image)

The extraneous lines are breaking up the interior land areas into multiple regions, which skews my total hole count to upwards of 270, instead of 221.

I researched all of the polyshape functions, looking for a way to clean up the data - such as "polybuffer" and "simplify". Apparently "simplify" runs by default when you do a merge. Perhaps my "exclusive or" is part of the problem? I tried zooming in on the line endpoints, looking for a Mandelbrot-like pocket or something, but even down to a thousandth of a degree resolution I see nothing. What do you suggest? It appears I am pushing the polyshape merge process beyond the normal limit.

[Kurt]

After further rumination... It appears that it makes more sense to compare the polyshape areas than to try counting holes. By this approach, comparing landmasses, the difference between my merged polygon and the "coast" polygon is only 1.6 square kilometers. I also got 22 holes, which is exactly what I should be seeing.

I also did NOT see any funny lines while aggregating landmasses. This suggests that perhaps the "xor" operation (only used for oceans) could be suspect.

Answer 1

[Matt J]

Perhaps I don't fully understand the logic of your loop completely, but I definitely think subtract() is more appropriate here than xor(). Also, if you are accumulating the oceans, it seems to me as though the union() should be taken with pgon_merge, not pgon.

for i = 1:numel(srcFiles)
    
    % read lat, lon from source files directory, 1 file at a time
    % lat, lon size: (286x636)
     
    
    filename = srcFiles(i).name;
    
    lat = ncread(filename,'Latitude');
    lon = ncread(filename,'Longitude');
    
    k = boundary(lat(:),lon(:));
    pgon = polyshape(lon(k),lat(k));
    
    
    
    pgon_sea =subtract(pgon,coast);		% sea only patch
    
    if i==1
       pgon_merge=pgon_sea;
    else
       pgon_merge = union(pgon_sea,pgon_merge);	% merge sea poly into global poly
    end
    
    plot(gx,pgon_merge);
    
    disp(strcat("patch holes = ",num2str(pgon_merge.NumHoles)));
    
end

[Matt J]

Note that a bug in union() was found almost a year ago that also caused region fractures,

https://www.mathworks.com/matlabcentral/answers/2115841-the-union-function-for-polyshapes-performs-an-incorrect-consolidation-of-adjacent-polyshapes-when?s_tid=srchtitle

Since you are using union() incrementally, I tend to doubt you are seeing the same bug, but it is something to be aware of all the same.

[Matt J]

I researched all of the polyshape functions, looking for a way to clean up the data - such as "polybuffer" and "simplify".

It does look to me like polybuffer could be necessary to make the union() robust if you have patches that closely border one another within whatever measurement error there is in your lat,lon data

tol=0.1; %tolerance against inexact lat,lon

pb=@(p,d) polybuffer(p,d);

pgon_merge = pb( union(  pb(pgon_sea,tol) ,  pgon_merge)   , -tol);

[Kurt]

Matt, thank you for these suggestions. I was not aware of the subtract() function. The pgon_merge thing is probably an error in my hand-copying of this code from another system.

I found a different approach where I compare the area of the polygons instead of counting holes. That way, I can ignore the fractured polygons - but measuring the area has its own challenges, especially around +/- 180 degrees longitude.

Some of the small polygons overlap at the ends, so as you say, polybuffer could be useful.

[Kurt]

Using "subtract" instead of "xor" appears to have eliminated the extraneous lines in my polyshapes. No other massaging with polybuffer was necessary. Thanks