Определить расстояние от береговой линии в Matlab

Question

В MATLAB у меня есть массив пар широты и долготы, которые представляют местоположения в Соединенных Штатах. Мне нужно определить расстояние до ближайшей береговой линии.

Я думаю, что в MATLAB есть встроенная база данных широт / долгот в Соединенных Штатах. Как я могу получить к нему доступ и использовать его?

Также есть предложения о том, как эффективно определить расстояние?

Обновление : следующий вопрос: Определить центр бункеров при использовании сетки

Answer 1

Поскольку у меня нет доступа к Mapping Toolbox , который идеально подходит для решения этой проблемы, я пришел к решению, которое не зависит от каких-либо наборов инструментов, включая Инструментарий обработки изображений .

Стив Эддинс имеет блог по обработке изображений на The MathWorks , где в прошлом году у него была серия довольно интересных постовпосвященный использованию цифровых карт рельефа.В частности, он указал, где их взять и как их загрузить и обработать.Вот соответствующие сообщения в блоге:

• Определение континентального водораздела США, часть 1 - Введение : Здесь Стив показывает, где можно получить цифровую карту рельефа (ЦМР)плитки и как их загружать и обрабатывать.Вы можете получить плитки DEM ( плитка E и плитка F охватывают континентальную часть США) из Global Land One-km Base Elevation Project .Диапазоны широты и долготы для каждой плитки можно найти здесь .

• Определение континентального водораздела США, часть 4 - Маски океана : Использованиеобработанные данные матрицы высот из вышеприведенного поста, Стив показывает, как создать маску океана.

С помощью этих данных матрицы высот вы можете узнать широту и долготу, где находятся края океанов.найдите расстояние от внутренней точки карты до ближайшей из этих прибрежных точек, а затем сделайте приятную визуализацию.Я использовал функцию FILTER2 , чтобы помочь найти края океанов посредством свертки с маской океана, и уравнения для вычисления расстояний большого круга , чтобы получить расстояние между точками карты вдольПоверхность Земли.

Используя некоторые примеры кода из вышеприведенных постов в блоге, я пришел к следующему:

%# 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 милях от ближайшей "океанской" береговой линии (на самом деле это, вероятно, Внутриреберный водный путь ).

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