MATLAB Answers

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

9 views (last 30 days)
Jamie LaPointe
Jamie LaPointe on 15 Feb 2019
Answered: Ted Shultz on 5 Jun 2020
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:
InverseAntarctica.png
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:
NYC.png
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.

Answers (2)

Jamie LaPointe
Jamie LaPointe on 9 May 2020
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

Ted Shultz
Ted Shultz on 5 Jun 2020
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);

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by