1

我对美国大陆许多点到海洋的 W、NW、SW 距离感兴趣。出于测试目的，我在 500 m（32x32 像素） GMTED2010和垂直海岸线上循环通过 1/8 度 dem 。我环顾了这个站点并因此实现了 pdist2 函数，但是我没有得到我期望的结果。所以我的第一个问题是我是否在概念上是错误的，第二个问题是我的 pdist2 实现是否不正确？我也对其他解决方案持开放态度。

考虑到方向约束，我希望在所有 3 个方向上看到相同的模式。最西边的像素列将具有相同的距离，下一列将是相同的，等等，所以当我绘制一个 32x32 矩阵dlong使用时，imagesc我得到一个从低到高、从左到右的渐变。

%**************
%For those truly interested, you can download the DEM and get Z and R accordingly:
[Z120,R120]=geotiffread('~/path/to/tif/GMTED2010N30W120_150/30n120w_20101117_gmted_mea150.tif');
[Z150,R150]=geotiffread('~/path/to/tif/GMTED2010N30W150_150/30n150w_20101117_gmted_mea150.tif');
Z=[Z150 Z120];
R=R120;
Z=Z(:,6001:4800+7200); %crop Z from -100 to -125. use latlon2pix to confirm between sub-z and z
R.Lonlim=[-125, -100]; 
R.RasterSize=size(Z);
clear Z150 Z120 R150 R120

%******* HERE STARTS THE ALGORITHM
%coastline (ultimately will be from the coast library)
latlim=[0.25:.25:60];
lonlim=ones(length(latlim),1)*-110

%variables r and c are the row and column indices for the point I'm interested in. r and c are relative to a DEM for the entire western USA so a point in Colorado is something like 2370,4350.
rstart=2370;
cstart=4350;

for r=2370:2370+31
    for c=4350:4350+31   
        %rows and cols are the vectors in the NW direction from point r,c.
        %in the SW direction, rows=r+[1:min(r,c)-1]. cols is the same.
        %W direction, rows=ones(r,1)*r; cols=c-[1:c-1];
        rows= r-[1:min(r,c)-1];
        cols= c-[1:min(r,c)-1];

        %Use referencing object R for DEM Z of the western USA to convert rows and cols to lat and long.
        [NWcoord(:,1) NWcoord(:,2)]=pix2latlon(R,rows,cols);

        %use pdist2 to find the shortest distance between any two points in the two vectors
        [D,i]=pdist2(lonlim,NWcoord(:,2),'euclidean','smallest',1);
        [~, mi]=min(D);

        sta.NWcoast=[latlim(i(mi)) lonlim(i(mi))];
        dlong(r-rstart+1,c-cstart+1)=distance(lat,long,latlim(i(mi)),lonlim(i(mi))); %great arc distance on earth's surface. radians
    end
end
4

2 回答 2

1

我的建议：

在没有DEM的情况下解决它。

1) 您已经给出了一个位置 loc，其纬度和经度以十进制度 WGS84 为单位。
2) 你已经在 WGS84 中给出了外套线多边形。

现在找到从 loc 到多边形的东西距离：
您想要找到 loc 的纬度值与边界多边形（当前位置以东）的交叉点：

从边界多边形点 0 的开头开始：找到border[i].lat<= loc.lat and border[i+1] > loc.lat AND border[i].longitude >= loc.longitude. 如果您找到了这条线，请在 (i 和 i+1) 之间进行线性插值，以找到准确的（纬度/经度）交点。

现在你有了与海洋的交点：计算距离 loc -> 与haversine公式的交点。

（一旦这工作，您可以稍后决定是否要加快二进制搜索）

与其他 3 个方向相同，交换 lat/long 和更大/更小

对于 NW 和其他： 沿着海洋边界点运行并计算从 loc 到边界点的方位（搜索航空公式或更大的圆方位计算）

存储线/或轴承跨过 315 度的两个点。然后这条线与 315° 相交，理论上可以有不止一条这样的线存储所有这样的线，并取离位置最近的一条（

现在用 315 对两个点进行插值以获得精确切割。

更新：轴承公式

于 2013-02-20T19:17:08.667 回答
0

如果您有参考对象 R，则以下工作有效。感谢@AlexWein 的方位角想法。之前发布的图片中出现条纹的原因是沿海矢量的比例。请注意，我使用的是 0.0042 度增量（与 dem 分辨率相同）。5 条条纹表示我只使用 5 个海岸坐标点来查找距离，而阶梯条纹的出现是因为我从下方最近点反弹到上方最近点并再次返回。

    latlim=[0.25:.0042:60];
    lonlim=ones(length(latlim),1)*-110

for r=1:32
    for c=1:32

    [lat long]=pix2latlon(R,r,c)
[d, az]=distance(lat,long,latlim,lonlim);

[~, azi]=min(abs(az-270));
sta_o.Wd2O=d(azi);
sta_o.Waz=az(azi);

[~, azi]=min(abs(az-315));
sta_o.NWd2O=d(azi);
sta_o.NWaz=az(azi);

[~, azi]=min(abs(az-225));
sta_o.SWd2O=d(azi);
sta_o.SWaz=az(azi);

    end
end
于 2013-03-12T00:55:14.827 回答