matlab - 在 Matlab 中确定距 coaSTLine 的距离

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 :在这里，史蒂夫展示了您可以获得数字高程 map (DEM) 切片的位置以及如何加载和处理它们。您可以从 tile E 获取 DEM 瓦片(tile F 和 Global Land One-km Base Elevation Project 覆盖美国大陆) .可以找到每个图 block 的纬度和经度范围 here .

• Locating the US continental divide, part 4 - Ocean masks :使用上述帖子中经过处理的 DEM 数据，史蒂夫展示了如何创建海洋掩模。

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

使用上述博客文章中的一些示例代码，这是我想出的:

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