投影计算

大多数情况下，人们只需要得到地图的图形表示，但有时候又确实需要地图数据在投影空间中的非地理坐标。[大谦MATLAB，dqmatlab点com]

用于提取地图投影坐标的简单方法之一就是使用MATLAB的get命令。投影坐标保存在对象的XData和YData属性中。

作为示例，下面显示Mollweide投影的图框，并提取出它的x-y坐标。

>> axesm mollweid
>> h=framem;
>> x=get(h,'XData');
>> y=get(h,'YData');

当然，不需要通过显示一个地图对象来获取它的投影坐标。可以用工具箱中的显示命令进行相同的投影计算。

下面从显示coast矢量数据的地图开始。

>> figure
>> load coast
>> plot(long,lat); axis equal
>> set(gca,'XLim',[-200 200],'YLim',[-100 100])

生成图4-11。

图4-11 coast矢量数据的地图

投影数据以前，必须定义投影参数。投影参数保存在一个地图投影结构中，该结构通常位于MATLAB axes对象的UserData属性中，但是可以直接用它进行投影计算。

下面的命令创建一个空的地图投影结构，改变地图原点，用结构的其他字段取默认值。

>> mstruct=defaultm('sinusoid');
>> mstruct.origin=[0 180 0];
>> mstruct=defaultm(sinusoid(mstruct));

定义地图投影参数以后，可以将经度矢量和纬度矢量投影到合适的坐标系中并用非映射的方法，即MATLAB命令获得结果。

>> [x,y]=mfwdtran(mstruct,lat,long,[],'line');
>> plot(x,y); axis equal
>> set(gca,'XLim',[-3.5 3.5],'YLim',[-2 2])

生成图4-12。

图4-12 重新定义投影参数以后的显示结果

在转换函数中，将纬度数据对应的参数设置为空矩阵，用'line'表示对象类型，并允许对矢量数据进行裁剪。

可以用逆转换函数将投影后的x-y数据转换回格林威治地理坐标。

>> [lat2,long2]=minvtran(mstruct,x,y);
>> plot(long2,lat2); axis equal
>> set(gca,'XLim',[-200 200],'YLim',[-100 100])

生成图4-13。

图4-13 将地图转换回格林威治坐标