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:
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).
¿Están esos archivos y funciones en la Caja de herramientas de asignación? – gnovice
Sí, estas son herramientas de mapeo. Estaba a punto de editar la respuesta para reflejar eso, ¡pero ya me atrapaste! – MatlabDoug
Tenía que hacer un diagrama (largo, lat) para ir al norte. – Marc