由于我无权访问绘图工具箱 http://www.mathworks.com/products/mapping/,这对于解决这个问题来说是理想的,我想出了一个独立于任何工具箱的解决方案,包括图像处理工具箱 http://www.mathworks.com/products/image/.
史蒂夫·埃丁斯 https://stackoverflow.com/users/172800/steve-eddins has an The MathWorks 的图像处理博客 http://blogs.mathworks.com/steve/去年,他发表了一系列非常酷的帖子,专门讨论如何使用数字高程地图。具体来说,他指出了从哪里获取它们以及如何加载和处理它们。以下是相关博客文章:
定位美国大陆分水岭,第 1 部分 - 简介 http://blogs.mathworks.com/steve/2009/04/24/continental-divide-1-intro/:在这里,Steve 展示了在哪里可以获得数字高程图 (DEM) 切片以及如何加载和处理它们。您可以获得 DEM 图块(tile E http://www.ngdc.noaa.gov/mgg/topo/DATATILES/elev/e10g.gz and tile F http://www.ngdc.noaa.gov/mgg/topo/DATATILES/elev/f10g.gz覆盖美国大陆)环球置地一公里基地高程项目 http://www.ngdc.noaa.gov/mgg/topo/globe.html。可以找到每个图块的纬度和经度范围here http://www.ngdc.noaa.gov/mgg/topo/gltiles.html.
定位美国大陆分水岭,第 4 部分 - 海洋面具 http://blogs.mathworks.com/steve/2009/05/15/continental-divide-4-oceans/:使用上述帖子中经过处理的 DEM 数据,Steve 展示了如何创建海洋掩模。
利用这些 DEM 数据,您可以找出海洋边缘所在位置的纬度和经度,找到从内陆地图点到最近的这些沿海点的距离,然后进行可视化。我使用了该功能FILTER2 http://www.mathworks.com/access/helpdesk/help/techdoc/ref/filter2.html通过与海洋掩模的卷积以及计算方程来帮助找到海洋的边缘大圆距离 http://en.wikipedia.org/wiki/Great-circle_distance获取沿地球表面的地图点之间的距离。
使用上述博客文章中的一些示例代码,我得出以下结论:
%# 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 英里(实际上,这可能是近岸水道 http://en.wikipedia.org/wiki/Intracoastal_Waterway).