Determine la distancia desde la costa en Matlab

Question

En MATLAB, tengo una serie de pares de latitud y longitud que representa ubicaciones en los Estados Unidos. Necesito determinar la distancia a la costa más cercana.

Creo que MATLAB tiene una base de datos integrada de los puntos lat/lon de los Estados Unidos. ¿Cómo accedo y lo uso?

¿Alguna sugerencia sobre cómo determinar de manera eficiente la distancia?

actualización: Seguimiento pregunta: Determine center of bins when using meshm

Answer 1

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

Hay otros conjuntos de datos en el mismo directorio que coast.mat que podría ser más útil.

Me gustaría encontrar la distancia a todos los puntos en el conjunto de datos, y tomar la distancia más corta. Esto supondría que las costas fuera de los EE. UU. Son respuestas aceptables. Querrá utilizar la función de distancia ya que la geometría euclidiana no se aplica aquí.

Answer 2

Ya que no tengo acceso a la Mapping Toolbox, lo que sería ideal para resolver este problema, se me ocurrió una solución que es independiente de cualquier cajas de herramientas, incluyendo el Image Processing Toolbox.

Steve Eddins tiene un image processing blog at The MathWorks donde el año pasado tuvo una serie de publicaciones muy interesantes dedicadas al uso de mapas de elevación digitales. Específicamente, señaló dónde obtenerlos y cómo cargarlos y procesarlos. Aquí están las entradas de blog relevantes:

• Locating the US continental divide, part 1 - Introduction: Aquí, Steve muestra donde se puede obtener un mapa de elevación digital azulejos (DEM) y cómo cargar y procesarlos. Puede obtener las fichas DEM (tile E y tile F en la cobertura continental de EE. UU.) Del Global Land One-km Base Elevation Project. Los rangos de latitud y longitud para cada mosaico se pueden encontrar here.

• Locating the US continental divide, part 4 - Ocean masks: Al usar los datos DEM procesados ​​de la publicación anterior, Steve muestra cómo crear una máscara oceánica.

Con estos datos DEM, se puede averiguar la latitud y longitud de donde los bordes de los océanos son, encontrar la distancia de un punto del mapa interior para el más cercano de estos puntos costeros, a continuación, hacer una visualización dulce . Usé la función FILTER2 para ayudar a encontrar los bordes de los océanos a través de la convolución con la máscara oceánica, y las ecuaciones para calcular great-circle distances para obtener la distancia entre los puntos del mapa a lo largo de la superficie de la Tierra.

El uso de algunos de los ejemplos de código de las entradas de blog por encima, esto es lo que ocurrió:

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

y aquí está la cifra resultante:

supongo que me pone casi 400 millas de la línea de costa "oceánica" más cercana (en realidad, es probablemente el Intracoastal Waterway).

Answer 3

La respuesta de Gnovice fue agradable y útil para el futuro, pero no necesitaba esa gran fidelidad y no quería gastar tiempo extra convirtiendo la distancia de píxel a la distancia de latitud/longitud. Tomando la respuesta de MatlabDoug, escribí la siguiente secuencia de comandos:

% 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