1. 多边形的分解和组装
当部分直线和面片被组合到元素用NaN分隔的很大的矢量中以后很难识别。可以用polysplit函数把这些多边形或直线矢量分解成基本图形单元,该函数将列矢量作为输入变量。[大谦MATLAB,dqmatlab点com]
下面的例子用polysplit函数提取和连接多边形或直线段。
输入两个列矢量形式的以NaN为界的数组。
>> lat=[45.6 -23.47 78 NaN 43.9 -67.14 90 -89]';
>> long=[13 -97.45 165 NaN 0 -114.2 -18 0]';
用polysplit函数创建两个单元数组latc和lonc。
>> [latc,lonc]=polysplit(lat,long)
latc=
[3x1 double] [4x1 double]
lonc=
[3x1 double] [4x1 double]
查看单元数组的内容。
>> [latc{1} lonc{1}]
ans =
45.6 13
-23.47 -97.45
78 165
>> [latc{2} lonc{2}]
ans =
43.9 0
-67.14 -114.2
90 -18
-89 0
注意,每个单元数组元素包含原直线的一段。
进行反操作,使用polyjoin函数。
>> [lat2,lon2]=polyjoin(latc,lonc);
连接后的线段与最初的lat和lon数组等价。
>> [lat long] == [lat2 lon2]
ans =
1 1
1 1
1 1
0 0
1 1
1 1
1 1
1 1
可以测试两套数据是否完全相等,键入下面的命令行。
>> isequalwithequalnans(lat,lat2) & isequalwithequalnans(long,lon2)
ans =
1
2. 将直线段连接成多边形
将具有公共端点的直线段进行连接是与直线段有关的一个通用操作。polymerge命令比较保存在纬度数据矢量和经度数据矢量中的线段端点数据,以便识别精确匹配或距离在一定范围内的端点,然后将匹配的直线段进行连接。重复此过程,直到没有公用点可以找到为止。有两个必需的变量是纬度矢量和经度矢量。
下面的例子通过将直线段连接成多边形演示了这个过程。
创建表示坐标值的列矢量。
>> lat=[3 2 NaN 1 2 NaN 5 6 NaN 3 4]';
>> lon=[13 12 NaN 11 12 NaN 15 16 NaN 13 14]';
将精确匹配的线段连接起来。
>> [latm,lonm]=polymerge(lat,lon)
ans =
5 15
6 16
NaN NaN
1 11
2 12
2 12
3 13
3 13
4 14
原来的4条线段变成了两条。
polymerge函数还有可选的第3个变量,它设置允许不精确匹配的距离容限。第4个变量使用户可以指定函数输出矢量还是单元数组。
3. 地理插值
可以用interpm函数在点之间进行线性填充。考虑一系列纬度和经度点,对它们进行线性插值处理,使得在每个方向上的距离都不超过1°。
>> lats=[1 2 4 5]; longs=[1 3 4 5]; maxdiff=1;
>> [newlats,newlongs]=interpm(lats,longs,maxdiff)
newlats =
1.0000
1.5000
2.0000
3.0000
4.0000
5.0000
newlongs =
1.0000
2.0000
3.0000
3.5000
4.0000
5.0000
在原来的lats数据中,2和4之间相差2,所以,newlats和newlongs两个新变量中插入了(3, 3.5)这个点。类似地,原来的longs数据中,1和3之间也相差2,所以newlats和newlongs两个变量中插入了(1.5, 2)这个点。现在,newlats和newlongs中再没有相邻的点的距离比maxdiff还大了。
用intrplat函数进行插值,可以选择用线性、样条、三次、恒向线或大圆弧等方法进行插值。利用下面给定的数据进行计算,求与经度7.3°对应的纬度,用线性、大圆和恒向线3种方法进行计算。
>> longs=[1 3 4 9 13]; lats=[57 68 60 65 56]; newlong=7.3;
>> newlat=intrplat(longs,lats,newlong,'linear')
newlat =
63.3000
>> newlat=intrplat(longs,lats,newlong,'gc')
newlat =
63.5029
>> newlat=intrplat(longs,lats,newlong,'rh')
newlat =
63.3937
intrplon函数的功能与intrplat的相似,只不过它求的是经度。
4. 矢量相交
Mapping工具箱提供了一系列函数对矢量数据进行求交计算。这些函数还可以计算任意矢量数据的交。
下面的例子计算圆心在(0°, 0°),半径为1250海里的小圆与圆心在(5ºN, 30ºE),半径为2500公里的小圆的交点。
>> [lat,long]=scxsc(0,0,nm2deg(1250),5,30,km2deg(2500))
lat =
17.7487 -12.9839
long =
11.0624 16.4170
其他功能相似的命令包括rhxrh、gcxgc和gcxsc等,其中,rhxrh函数计算恒向线的交点,gcxgc函数计算大圆的交点,gcxsc函数计算大圆和小圆的交点。
5. 多边形的面积
可以用函数areaint计算多边形格式矢量数据的地理面积。
下面用usalo数据集计算美国大陆的面积。
>> load usalo
>> earthradius=almanac('earth','radius');
>> area=areaint(uslat,uslon,earthradius)
area =
1.0e+06 *
7.9256
0.0035
0.0004