Plotting Antarctica filled polygon with gshhs function fills the inverse of Antarctica

I am using the gshhs function to load the latest version of GSHHG coastal data (version 2.3.7). This data was downloaded from https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/. Information about this dataset: https://www.ngdc.noaa.gov/mgg/shorelines/

I plot the filled polygons using geoshow and for the parts of the Earth north of S (e.g. Anarctica) it works great. However, for Antarctica, it fills the inverse of Anartica instead.

Here is example code use to plot Antarctica:

% Antarctica test

figure('Name', 'Antartica', 'NumberTitle', 'off', 'Color', 'white');
worldmap('Antarctica');

projection = gcm;
latlim = projection.maplatlimit;
lonlim = projection.maplonlimit;

if ~exist('gshhs_f.i', 'file');
    gshhs('gshhs_f.b', 'createindex');
end

% Load the GSHHG coastal polygon data version 2.3.7 - Full Resolution
antarctica = gshhs('gshhs_f.b', latlim, lonlim);

levels          = [antarctica.Level];
land            = (levels == 1);
lake            = (levels == 2);
island          = (levels == 3); % island in a lake
pond            = (levels == 4); % pond in an island in a lake
ice_front       = (levels == 5); % ice shelves around Antarctica
grounding_line  = (levels == 6); % land of Antarctica

geoshow([antarctica(ice_front).Lat],      [antarctica(ice_front).Lon],      'DisplayType', 'Polygon', 'FaceColor', [230/255 230/255 230/255]); % gray
geoshow([antarctica(grounding_line).Lat], [antarctica(grounding_line).Lon], 'DisplayType', 'Line',    'Color',     [255/255 105/255 180/255]); % hot pink
geoshow([antarctica(land).Lat],           [antarctica(land).Lon],           'DisplayType', 'Polygon', 'FaceColor', [  0/255 100/255   0/255]); % forest green
geoshow([antarctica(lake).Lat],           [antarctica(lake).Lon],           'DisplayType', 'Polygon', 'FaceColor', [  0/255   0/255 128/255]); % navy blue
geoshow([antarctica(island).Lat],         [antarctica(island).Lon],         'DisplayType', 'Polygon', 'FaceColor', [210/255 105/255  30/255]); % chocolate
geoshow([antarctica(pond).Lat],           [antarctica(pond).Lon],           'DisplayType', 'Polygon', 'FaceColor', [ 84/255  84/255  84/255]); % light steel blue

setm(gca, 'FFaceColor', 'cyan'); % set background color (should be seas, oceans, bays, etc.)

box('off');
xlabel('lon (deg)', 'fontsize', 12, 'color', 'b');
ylabel('lat (deg)', 'fontsize', 12, 'color', 'b');
axis('tight');
axis('fill');
grid('on');

The results of said plot:

The light gray is supposed to be Antarctica and the cyan is supposed to be the seas, oceans, bays, etc. Clearly it is the inverse.

Here is a simple example of another part of the world (New York City, NY, USA) that plots just fine:

% NYC test

figure('Name', 'NYC, USA', 'NumberTitle', 'off', 'Color', 'white');

latlim = [40 41];
lonlim = [-75 -73];
worldmap(latlim, lonlim);

% Load the GSHHG coastal polygon data version 2.3.7 - Full Resolution
if ~exist('gshhs_f.i', 'file');
    gshhs('gshhs_f.b', 'createindex');
end
nyc = gshhs('gshhs_f.b', latlim, lonlim);

levels          = [nyc.Level];
land            = (levels == 1);
lake            = (levels == 2);
island          = (levels == 3); % island in a lake
pond            = (levels == 4); % pond in an island in a lake
ice_front       = (levels == 5); % ice shelves around Antarctica
grounding_line  = (levels == 6); % land of Antarctica

geoshow([nyc(ice_front).Lat],      [nyc(ice_front).Lon],      'DisplayType', 'Polygon', 'FaceColor', [230/255 230/255 230/255]); % gray
geoshow([nyc(grounding_line).Lat], [nyc(grounding_line).Lon], 'DisplayType', 'Line',    'Color',     [255/255 105/255 180/255]); % hot pink
geoshow([nyc(land).Lat],           [nyc(land).Lon],           'DisplayType', 'Polygon', 'FaceColor', [  0/255 100/255   0/255]); % forest green
geoshow([nyc(lake).Lat],           [nyc(lake).Lon],           'DisplayType', 'Polygon', 'FaceColor', [  0/255   0/255 128/255]); % navy blue
geoshow([nyc(island).Lat],         [nyc(island).Lon],         'DisplayType', 'Polygon', 'FaceColor', [210/255 105/255  30/255]); % chocolate
geoshow([nyc(pond).Lat],           [nyc(pond).Lon],           'DisplayType', 'Polygon', 'FaceColor', [ 84/255  84/255  84/255]); % light steel blue

setm(gca, 'FFaceColor', 'cyan'); % set background color (should be seas, oceans, bays, etc.)

box('off');
xlabel('lon (deg)', 'fontsize', 12, 'color', 'b');
ylabel('lat (deg)', 'fontsize', 12, 'color', 'b');
axis('tight');
axis('fill');
grid('on');

Here is its output:

This was the output expected.

I understand this is probably an issue related to the singularites near the poles and MATLAB not quite understanding how to draw a polygon that is open in rectilinear space but closed in spherical space. I have attempted some simple attempts at closing the polygon in rectilinear space myself to no avail. Any assistance with plotting this data in MATLAB R2017b using the full resolution of the GSSHG coastal data set would be much appreciated.

Answer 1

[Jamie LaPointe]

Use the knowledge of the reversal to arrive at a simple solution

close all; clear all; %#ok

constants = terrain_constants();

latlim = [-90 -50];
lonlim = [-180 180];

shorelines = gshhs(constants.gshhg_file, latlim, lonlim); % full-resolution GSHHG 2.3.7 binary file

shoreline_land       = shorelines([shorelines.Level]==constants.gsshg_level_land); % land == 1
shoreline_antarctica = shorelines([shorelines.Level]==constants.gsshg_level_ice_front); % ice_front == 5

% Earth sans Antarctica
lat = [shoreline_land.Lat];
lon = [shoreline_land.Lon];

% Antarctica
lat1 = [shoreline_antarctica.Lat];
lon1 = [shoreline_antarctica.Lon];

% Earth sans Antarctica:
%  0 == land
%  1 == coast line
%  2 == ocean/sea
[Z0,R0] = vec2mtx(lat, lon, 120, latlim, lonlim, 'filled');

% Antarctica is the reverse...:
%  0 == ocean/sea
%  1 == coast line
%  2 == land
[Z1,R1] = vec2mtx(lat1, lon1, 120, latlim, lonlim, 'filled');

Z = Z0; R = R0;
Z(Z0 == 2 & Z1 == 2) = 3; % or set Antarctica to 1 to just make it land
% Z now is:
%   0 == land sans Antarctica
%   1 == coast line
%   2 == ocean/sea
%   3 == Antarctica

% Plot the coastline with filled in colours
figure('Name',        'GSHHG Coastline', ...
       'NumberTitle', 'off',             ...
       'MenuBar',     'none',            ...
       'Color',       'white');
worldmap(latlim,lonlim);
geoshow(Z,R,'DisplayType','texturemap');
colormap([0 100/255 0;              ... % land sans Antarctica = forest green
          0 0 0;                    ... % coastline  = black
          0 1 1;                    ... % ocean/sea  = cyan
          255/255 105/255 180/255]);    % Antarctica = hot pink

Answer 2

[Ted Shultz]

I had the exact same problem, but found a simple solution. Flip the order of the points in Antarctica!

The problem is that geoshow assumes that the data is clockwise (or counter clockwise, I cant remember which), and in the GSHHG data set, it is the opposite for the continent of Antarctica.

I think in my dataset, the continent of Antarctica was the 5 object, so my hacky fix was:

antarcticaDat = dataSet(5);
antarcticaDat = fliplr(antarcticaDat.Lat);
antarcticaDat = fliplr(antarcticaDat.Lon);