在 Matlab 中确定与海岸线的距离

Question

在 MATLAB 中，我有一组表示美国位置的纬度和经度对。我需要确定到最近海岸线的距离。

我认为 MATLAB 有一个内置的美国纬度/经度点数据库。如何访问和使用它？

还有关于如何有效地确定距离的任何建议？

更新：后续问题：Determine center of bins when using meshm

Answer 1

由于我无法访问 Mapping Toolbox，这对于解决这个问题来说是理想的，我想出了一个独立于任何工具箱的解决方案，包括 Image Processing Toolbox。

Steve Eddins 有一个image processing blog at The MathWorks，去年他在其中发表了一系列非常酷的帖子，专门讨论使用数字高程图。具体来说，他指出了从哪里获取它们以及如何加载和处理它们。以下是相关的博文：

• Locating the US continental divide, part 1 - Introduction：在这里，Steve 展示了在哪里可以获得数字高程图 (DEM) 切片以及如何加载和处理它们。您可以从Global Land One-km Base Elevation Project 获取 DEM 瓦片（tile E 和 tile F 覆盖美国大陆）。每个图块的经纬度范围可以在here找到。

• Locating the US continental divide, part 4 - Ocean masks：使用上述帖子中处理后的 DEM 数据，Steve 展示了如何创建海洋掩膜。

使用此 DEM 数据，您可以找出海洋边缘所在位置的纬度和经度，找出从内陆地图点到最近这些沿海点的距离，然后进行精美的可视化。我使用函数FILTER2通过与海洋掩码的卷积来帮助找到海洋的边缘，并使用计算方程great-circle distances来获得沿地球表面的地图点之间的距离。

使用上述博客文章中的一些示例代码，我得出以下结论：

%# 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

这是结果图：

我想这让我距离最近的“海洋”海岸线近 400 英里（实际上，它可能是Intracoastal Waterway）。

Answer 2

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

在与coast.mat 相同的目录中还有其他数据集可能更有用。

然后我会找到到数据集中所有点的距离，然后取最短距离。这将假设美国以外的海岸线是可以接受的答案。您将需要使用距离函数，因为欧几里得几何在这里不适用。

Answer 3

Gnovice 的回答很好，对未来很有用，但我不需要那么高的保真度，也不想花费额外的时间从像素距离转换为纬度/经度距离。以 MatlabDoug 的回答，我编写了以下脚本：

% Get Data  
coast = load('coast.mat');
locations = load('locations.mat');

% Preallocate  
coast_indexes = nan(size(locations.lat));
distancefromcoast = nan(size(locations.lat));

% Find distance and corresponding coastal point  
for i=1:1:numel(locations.lat)  
    [dist, az] = distance(locations.lat(i), locations.long(i), coast.lat, coast.long);
    [distancefromcoast(i),coast_indexes(i)] = min(dist);
end