matlab example
Written by Daniel Herber on June 12, 2017.
While I was exploring [Google Maps](https://www.google.com/maps/@26.3444557,-40.9843087,8277695m/data=!3m1!1e3), I started to wonder what point in the ocean was the farthest away from land. So I fired up MATLAB to determine this location and visualize the distance to the nearest coast for each point in the ocean.
Later in this process, I found the proper name for this point: [oceanic pole of inaccessibility](https://en.wikipedia.org/wiki/Pole_of_inaccessibility). According to the Wikipedia page, it is located at [(48°52.6′S 123°23.6′W)](https://tools.wmflabs.org/geohack/geohack.php?pagename=Pole_of_inaccessibility¶ms=48_52.6_S_123_23.6_W_&title=Oceanic+Pole+of+Inaccessibility).
---
### Result
[![map_opi.png](blogs/matlab/post_7/map_opi.png)](blogs/matlab/post_7/map_opi.png){data-lightbox="blog_imgs" data-title="Oceanic Pole of Inaccessibility"}
Using a 1600 x 800 grid of points and the high-resolution dataset, the oceanic pole of inaccessibility was found to be (48°59.9′S 123°29.4′W) which is an error of about 9 miles from the true location.
---
### The Code
The code is broken down into three steps: 1) obtaining the data; 2) calculating the distances; 3) plotting the results. There are also three user-parameters:
```matlab
level = 4; % detail level of the map, between 1-5 (see obtaining the data)
N = 800; % number of latitude test points (see calculating the distances)
Ncontours = 256; % number of contours (see plotting the results)
% ensure that the data and export_fig are in your path
```
The complete code can be downloaded from [[link]](blogs/matlab/post_7/visualizing_opi.m).
(Some additional post-processing was done in [Inkscape](https://inkscape.org) to add borders)
#### Obtaining the Data
A good source for coastline data is the [GSHHG](https://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html) or Global Self-consistent, Hierarchical, High-resolution Geography Database.
For this project, I downloaded [gshhg-shp-2.3.6.zip](https://www.ngdc.noaa.gov/mgg/shorelines/data/gshhg/latest/) which contains ```.shp``` files. This data comes in five resolutions from crude to full. Lower resolutions were used when rapid prototyping the code and the highest resolution used was high since the full data proved challenging to work with. Antarctica is handled a bit differently than the other continents since the coastline can either be defined by the ice front or grounding line boundary. Other work calculating the oceanic pole of inaccessibility seemed to use the grounding line boundary, so it was selected for this project as well.
With the ```.shp``` files in our path, the data is then loaded with:
```matlab
switch level
case 1 % Crude resolution
levelstr = 'c'; % 8,526 points
case 2 % Low resolution
levelstr = 'l'; % 64,870 points
case 3 % Intermediate resolution
levelstr = 'i'; % 385,335 points
case 4 % High resolution
levelstr = 'h'; % 1,820,922 points
case 5 % Full resolution
levelstr = 'f'; % 9,985,362 points
end
% load continental land masses and ocean islands, except Antarctica
S1 = shaperead(['GSHHS_',levelstr,'_L1.shp'], 'UseGeoCoords', true);
% load Antarctica based on grounding line boundary
S6 = shaperead(['GSHHS_',levelstr,'_L6.shp'], 'UseGeoCoords', true);
% combine
coastLAT = [ [S1.Lat], [S6.Lat] ]; % coast latitude data
coastLON = [ [S1.Lon], [S6.Lon] ]; % coast longitude data
```
#### Calculating the Distances
The next step was to determine the test points and calculate the distance to the nearest coastline. First, we find the test points that are on land using the [inpolygon](https://www.mathworks.com/help/matlab/ref/inpolygon.html) for each of the polygons in the dataset. Next, for each test point that is not on land, we calculate the great circle arcs connecting pairs of points on the surface of a sphere using [distance](https://www.mathworks.com/help/map/ref/distance.html). Finally, we convert from degrees to miles using [distdim](https://www.mathworks.com/help/map/ref/distdim.html). While the distances calculated, the maximum value can be found to determine the oceanic pole of inaccessibility.
The code for this step is:
```matlab
% bounds on lat and long (degrees)
minLON = -180; maxLON = 180; minLAT = -90; maxLAT = 90;
% create grid of test points
LONv = linspace(minLON, maxLON, 2*N); % equally spaced longitude points
LATv = linspace(minLAT, maxLAT, N); % equally spaced latitude points
[testLON,testLAT] = meshgrid(LONv, LATv);
% determine the points that are on land
keep = ones(N, 2*N); % initialize as not inside a land mass
polyI = [0,find(isnan(coastLAT))]; % find the boundaries of the polygons
for idx = 1:length(polyI)-1
localI = polyI(idx)+1:polyI(idx+1)-1; % indices of the current polygon
% determine if the test points are inside the current polygon
in = inpolygon(testLAT(:), testLON(:), coastLAT(localI), coastLON(localI));
keep(in) = 0; % zero points inside the current polygon
end
% initialize and calculate distances
DIST = zeros(N,2*N);
parfor idx = 1:2*N^2
if keep(idx) % only if the point is in water
DIST(idx) = min(distance(testLAT(idx), testLON(idx), coastLAT, coastLON));
end
end
% convert degrees to miles
DISTmi = distdim(DIST,'deg','miles');
% find the oceanic pole of inaccessibility
[~,minI] = max(DIST(:));
```
#### Plotting the Results
The final step is to visualize the results using the [mapping toolbox](https://www.mathworks.com/products/mapping.html). A worldmap is initialized using [worldmap](https://www.mathworks.com/help/map/ref/worldmap.html). Next, the contours for the distances are drawn using [contourfm](https://www.mathworks.com/help/map/ref/contourfm.html). The land masses are plotted in black using [geoshow](https://www.mathworks.com/help/map/ref/geoshow.html) with the polygon display type. The oceanic pole of inaccessibility can then be plotted and a circle with radius equal to the distance to the nearest coastline is drawn utilizing [scircle1](https://www.mathworks.com/help/map/ref/scircle1.html) and [plotm](https://www.mathworks.com/help/map/ref/plotm.html). A colorbar and title are then added. Finally, the figure is saved as a ```.png``` using the [export_fig](https://github.com/altmany/export_fig) package.
The code for this step is:
```matlab
hf = figure; % create a new figure and save handle
hf.Color = [1 1 1]; % change the figure background color
hf.Position = [0 0 1920 1080]; % set figure size and position
ax = worldmap('World'); % initialize world map
gridm('off'); mlabel('off'); plabel('off'); % turn off labels and lines
ax.Children(end).EdgeColor = [0.4 0.4 0.4]; % change border color
hold on % hold the current graph
% plot the distance contours
colormap(flipud(winter(Ncontours))) % change the colormap
[cs,hc] = contourfm(testLAT, testLON, DISTmi,...
linspace(0,max(DISTmi(:)),Ncontours), 'LineColor', 'none');
% plot land (black)
geoshow([S1.Lat], [S1.Lon], 'DisplayType', 'polygon', 'FaceColor', [0 0 0]);
geoshow([S6.Lat], [S6.Lon], 'DisplayType', 'polygon', 'FaceColor', [0 0 0]);
% plot the oceanic pole of inaccessibility
redcolor = [244 67 54]/255; % material red color
geoshow(testLAT(minI), testLON(minI), 'DisplayType', 'point',...
'markersize', 14, 'MarkerEdgeColor', redcolor, 'linewidth', 2); % marker
[latc,longc] = scircle1(testLAT(minI), testLON(minI), DIST(minI));
plotm(latc, longc, 'color', redcolor, 'linewidth', 2); % circle
% add a colorbar
hcolorbar = colorbar;
hcolorbar.FontName = 'Segoe UI'; hcolorbar.FontSize = 18;
hcolorbar.TickLabels = strcat(hcolorbar.TickLabels,' mi'); % mile label
hcolorbar.Color = [0 0 0]; % black
% add title
htitle = title('Oceanic Pole of Inaccessibility');
htitle.FontName = 'Segoe UI';
htitle.FontSize = 32; htitle.FontWeight = 'normal';
% save the figure as a png
str = ['export_fig ''','map_opi',''' -png -m1 -transparent -opengl'];
eval(str)
```