How to prevent wraparound using geoplot

I am trying to accomplish something similar to the below post, namely prevent my data from wrapping at the meridians. So far, no luck:

https://www.mathworks.com/matlabcentral/answers/2172159-how-to-wrap-axes-using-geoplot

If I set the lonLim to [-360 360] and call geolimits(latlim,lonlim), the longitude limits do NOT change to [-360 360] but to [-180.3914 181.1743]. Why is this? Is there a bult-in max limit in geolimits? It's possible that the function is confused by the [-360 360] values, which are effectively equivalent.

But back to the wraparound issue. I am plotting NetCDF-type satellite observation data.

figure
gx = geoaxes;
geobasemap(gx,'usgsimagery');
hold on
latLim = [-72 85];
lonLim = [-180 180];
geolimits(gx, latLim, lonLim);

lat = ncread(filepath, 'Latitude');     % read NetCDF-type data
lon = ncread(filepath, 'Longitude');
lat1 = double(lat(:,:,1)) ;             % read just the 1st of 6 layers
lon1 = double(lon(:,:,1));
lon1 = wrapTo180(lon1);                 % limit to +/- 180 degrees
k = boundary(lat1(:),lon1(:));          % get just the edge data
pgon = geopolyshape(lat1(k),lon1(k));
geoplot(gx,pgon);

I am able to plot the satellite observation patches on the Mercator projection geoaxes display. However, when a patch crosses the International Date Line, it does not straddle the line but whangs all the way to the right, to the Prime Meridian.

If I use geoplot instead of geopolyshape I get the same effect, only with vectors instead of shapes.When the longitude data reaches -180, the vectors jump all the way to the right instead of straddling the dateline.

geoplot(gx,lat1(k),lon1(k));

How can I massage my NetCDF data to prevent this?

My goal is to collect all the patches and merge them into a single surface as described here:

https://www.mathworks.com/matlabcentral/answers/80524-combine-surfaces-into-one-big-surface

[Kurt]

The longitude tables in question have a discontinuity in the column data:

-147.75 -152.944 -158.91 -165.725 -173.292 178.5214 170.0151 161.5731 153.4542 146.0499

Somehow, I need to rectify this.

[Walter Roberson]

Deg = [-147.75 -152.944 -158.91 -165.725 -173.292 178.5214 170.0151 161.5731 153.4542 146.0499 ]

Deg = 1×10

-147.7500 -152.9440 -158.9100 -165.7250 -173.2920 178.5214 170.0151 161.5731 153.4542 146.0499

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

rad2deg(unwrap(deg2rad(Deg)))

ans = 1×10

-147.7500 -152.9440 -158.9100 -165.7250 -173.2920 -181.4786 -189.9849 -198.4269 -206.5458 -213.9501

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Walter Roberson]

Deg = [-147.75 -152.944 -158.91 -165.725 -173.292 178.5214 170.0151 161.5731 153.4542 146.0499 ]

Deg = 1×10

-147.7500 -152.9440 -158.9100 -165.7250 -173.2920 178.5214 170.0151 161.5731 153.4542 146.0499

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Deg(Deg<0) = 360 + Deg(Deg<0)

Deg = 1×10

212.2500 207.0560 201.0900 194.2750 186.7080 178.5214 170.0151 161.5731 153.4542 146.0499

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Kurt]

It's working better after I added the unwrap function, but I still have issues around the Prime Meridian. Patch data to the east of Greenwich looks OK. Data to the west is mirror-imaged and is still whanging 360 degrees east across the map to the International Date Line, instead of showing up normally to the west. The lonLim is [-180 180] but this doesn't appear to affect what I'm doing.

I failed to point out earlier that I am working with 2d polygons, not just vectors. lat1 and lon1 are 571x318 double arrays. When I used unwrap() I specified it as follows:

latUnwrap = rad2deg(unwrap(deg2rad(lat1),[],2));
lonUnwrap = rad2deg(unwrap(deg2rad(lon1),[],2));
k = boundary(latunwrap(:), lonunwrap(:));
pgon = geopolyshape(latunwrap(k),lonunwrap(k));
geoplot(gx,pgon);

I assume I have to unwrap latitude as well, so the lat/lon pairs in the two arrays still match up?

If all else fails, I suppose I can look for polygons that straddle the Prime Meridian and split them into two polygons, one on each side of the line?

This is purely a display problem, and does not affect the polygons themselves.

[Kurt]

(Later) I found the problem, but I'm not sure yet what the solution is. This is another discontinuity, but in the vertical direction in the longitude table:

359.9445 359.9153 359.8865 ...
359.9788 359.9496 359.9208 ...
0.013182 359.9839 359.9551 ...
0.047721 0.018456 359.9896 ...

This is right along the Prime Meridian, and instead of plotting to the left of the meridian it draws in the other direction, clear across the view to the right opposite edge. Solution? Another unwrap, perhaps?

Answer 1

[Kurt]

The solution I came up with was to split the patch in two along the meridian: half to the left, and half to the right. There may be cleaner and simpler approaches, but this worked.

latLeft = [];
lonLeft = [];
latRight = [];
lonRight = [];
m = [];
n = [];
left = false;
right = false;
lonUnwrap180 = mod((lonUnwrap+180),360) - 180;

for j = 1:height(lonUnwrap)
    idx = any(lonUnwrap(j,:) < 0);  % west of prime meridian
    idy = any(lonUnwrap(j,:) > 0);  % east of prime meridian
    if idx == 1
        lonLeft = [lonLeft;lonUnwrap(j,:)];
        m = [m,idx];        % build index to trim lat left
        left = true;
    end
    if idy == 1
        lonRight = [lonRight;lonUnwrap(j,:)];
        n = [n;idy];        % build index to trim lat right
        right = true;
    end
end

if left
    latLeft = latUnwrap(m==1,:);    % trim latitude to match
else
    latLeft = [];           % nothing to split
    lonLeft = [];   
end
if right
    latRight = latUnwrap(n==1,:)
    if ~left
        lonRight = lonUnwrap;
    else
        lonRight = lonUnwrap180;     % split - left & right
    end
else
    latRight = latUnwrap;
end

if left
    k = boundary(latLeft(:),lonLeft(:));    % get the edge data
    pgon = geopolyshape(latLeft(k),lonLeft(k));
    geoplot(gx,pgon);
end

k = boundary(latRight(:),lonRight(:));  % right is default if no split
pgon = geopolyshape(latRight(k),lonRight(k));
geoplot(gx,pgon);