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

5 次查看（过去 30 天）
Jamie LaPointe
Jamie LaPointe 2019-2-15
回答： Ted Shultz 2020-6-5
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.

请先登录，再回答此问题。

回答（2 个）

Jamie LaPointe
Jamie LaPointe 2020-5-9
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 2020-6-5
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);

类别

在 Help Center 和 File Exchange 中查找有关 Map Customization 的更多信息

标签

产品


版本

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by