1

米国本土の多くの地点の海までの西、北西、南西の距離に興味があります。テスト目的で、500 m (32x32 ピクセル) GMTED2010と垂直海岸線で 1/8 度の dem をループしています。このサイトを見回した結果、 pdist2 関数を実装しましたが、期待どおりの結果が得られません。だから私の最初の質問は、私が概念的に間違っているかどうかであり、2 番目の質問は私の pdist2 実装が間違っているかどうかです。私は他の解決策にもオープンです。

方向の制約が与えられた場合、3 方向すべてで同じパターンが見られると思います。ピクセルの最も西の列は同じ距離になり、次の列は同じになります。そのため、dlong使用する 32x32 マトリックスをプロットする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）10進角WGS84で緯度と経度を使用してロケーションlocを指定しました。
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方向についても同じように、緯​​度/経度と大/小を交換します。

NWおよびその他の場合： 海の境界点に沿って実行し、位置から境界点までの方位を計算します（航空の公式またはより大きな円方位の計算を検索します）

方位が315度を超えるライン/または2つのポイントを格納します。次に、この線は315°と交差します。そのような線は複数存在する可能性があり、そのような線をすべて格納し、場所に最も近い1つの線を取得します（

次に、両方のポイントを補間して、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 に答える