Determine distance from coastline in Matlab

Question

In MATLAB, I have an array of latitude and longitude pairs that represents locations in the United States. I need to determine the distance to the nearest coastline.

I think MATLAB has a built in database of lat/lon points of the United States. How do I access it and use it?

Also any suggestions on how to efficiently determine the distance?

Update: Follow-up question : Determine center of bins when using meshm

Answer 1

Since I don't have access to the Mapping Toolbox, which would be ideal to solve this problem, I came up with a solution that is independent of any toolboxes, including the Image Processing Toolbox.

Steve Eddins has an image processing blog at The MathWorks where last year he had a series of pretty cool posts devoted to using digital elevation maps. Specifically, he pointed out where to get them and how to load and process them. Here are the relevant blog posts:

• Locating the US continental divide, part 1 - Introduction: Here, Steve shows where you can get digital elevation map (DEM) tiles and how to load and process them. You can get the DEM tiles (tile E and tile F cover the continental US) from the Global Land One-km Base Elevation Project. The latitude and longitude ranges for each tile can be found here.

• Locating the US continental divide, part 4 - Ocean masks: Using the processed DEM data from the above post, Steve shows how to create an ocean mask.

With this DEM data, you can find out the latitude and longitude for where the edges of the oceans are, find the distance from an inland map point to the closest of these coastal points, then make a sweet visualization. I used the function FILTER2 to help find the edges of the oceans via convolution with the ocean mask, and the equations for computing great-circle distances to get the distance between map points along the surface of the Earth.

Using some of the sample code from the above blog posts, here's what I came up with:

%# Load the DEM data:

data_size = [6000 10800 1];  %# The data has 1 band.
precision = 'int16=>int16';  %# Read 16-bit signed integers into a int16 array.
header_bytes = 0;
interleave = 'bsq';          %# Band sequential. Not critical for 1 band.
byte_order = 'ieee-le';
E = multibandread('e10g',data_size,precision,...        %# Load tile E
                  header_bytes,interleave,byte_order);
F = multibandread('f10g',data_size,precision,...        %# Load tile F
                  header_bytes,interleave,byte_order);
dem = [E F];  %# The digital elevation map for tile E and F
clear E F;    %# Clear E and F (they are huge!)

%# Crop the DEM data and get the ranges of latitudes and longitudes:

[r,c] = size(dem);      %# Size of DEM
rIndex = [1 4000];      %# Row range of DEM to keep
cIndex = [6000 14500];  %# Column range of DEM to keep
dem = dem(rIndex(1):rIndex(2),cIndex(1):cIndex(2));  %# Crop the DEM
latRange = (50/r).*(r-rIndex+0.5);     %# Range of pixel center latitudes
longRange = (-180/c).*(c-cIndex+0.5);  %# Range of pixel center longitudes

%# Find the edge points of the ocean:

ocean_mask = dem == -500;        %# The ocean is labeled as -500 on the DEM
kernel = [0 1 0; 1 1 1; 0 1 0];  %# Convolution kernel
[latIndex,longIndex] = ...       %# Find indices of points on ocean edge
  find(filter2(kernel,~ocean_mask) & ocean_mask);
coastLat = latRange(1)+diff(latRange).*...     %# Convert indices to
           (latIndex-1)./diff(rIndex);         %#   latitude values
coastLong = longRange(1)+diff(longRange).*...  %# Convert indices to
            (longIndex-1)./diff(cIndex);       %#   longitude values

%# Find the distance to the nearest coastline for a set of map points:

lat = [39.1407 35 45];        %# Inland latitude points (in degrees)
long = [-84.5012 -100 -110];  %# Inland longitude points (in degrees)
nPoints = numel(lat);         %# Number of map points
scale = pi/180;               %# Scale to convert degrees to radians
radiusEarth = 3958.76;        %# Average radius of Earth, in miles
distanceToCoast = zeros(1,nPoints);   %# Preallocate distance measure
coastIndex = zeros(1,nPoints);        %# Preallocate a coastal point index
for iPoint = 1:nPoints                %# Loop over map points
  rho = cos(scale.*lat(iPoint)).*...  %# Compute central angles from map
        cos(scale.*coastLat).*...     %#   point to all coastal points
        cos(scale.*(coastLong-long(iPoint)))+...
        sin(scale.*lat(iPoint)).*...
        sin(scale.*coastLat);
  d = radiusEarth.*acos(rho);         %# Compute great-circle distances
  [distanceToCoast(iPoint),coastIndex(iPoint)] = min(d);  %# Find minimum
end

%# Visualize the data:

image(longRange,latRange,dem,'CDataMapping','scaled');  %# Display the DEM
set(gca,'DataAspectRatio',[1 1 1],'YDir','normal',...   %# Modify some axes
    'XLim',longRange,'YLim',fliplr(latRange));          %#   properties
colormap([0 0.8 0.8; hot]);  %# Add a cyan color to the "hot" colormap
xlabel('Longitude');         %# Label the x axis
ylabel('Latitude');          %# Label the y axis
hold on;                     %# Add to the plot
plot([long; coastLong(coastIndex).'],...    %'# Plot the inland points and
     [lat; coastLat(coastIndex).'],...      %'#   nearest coastal points
     'wo-');
str = strcat(num2str(distanceToCoast.',...  %'# Make text for the distances
                     '%0.1f'),{' miles'});
text(long,lat,str,'Color','w','VerticalAlignment','bottom');  %# Plot the text

And here's the resulting figure:

I guess that puts me almost 400 miles from the nearest "ocean" coastline (in actuality, it's probably the Intracoastal Waterway).

Answer 2

load coast; 
axesm('mercator'); 
plotm(lat,long)

There are other datasets in the same directory as coast.mat that might be more useful.

I would then just find the distance to all the points in the dataset, and take the shortest distance. This would assume that coastlines outside of the US are acceptable answers. You will want to use the distance function since Euclidean geometry does not apply here.