摘 要：在MATLAB平台下，综合利用等差分纬线多圆锥投影的正反解变换实现了投影后矢量和栅格数据的输出。输出后的成图不仅保证了数据的准确性，可靠性，并且可以供用户直接在地理信息系统软件ArcMap中读取和编辑，方便用户直接使用。
关键词：等差分纬线多圆锥投影 MATLAB 投影变换 TFW
中图分类号：P226 文献标识码：A 文章编号：1672-3791(2013)06(b)-0039-03
等差分纬线多圆锥投影是一种属于任意性质的多圆锥投影，由我国地图工作者于1963年根据我国同周围的国家和地区的位置和形状关系自主设计的一种投影。使用这种投影制作世界地图时，往往选取东经150°作为中央经线，这样不仅能完整显示太平洋和各洲大陆，还能使我国版图处于图上的相对居中位置[1，2]。因此该投影在中国大陆出版的各种世界地图专题图中被广泛应用，取得了较好的效果。
等差分纬线多圆锥投影的性质为：投影后的纬线(包括极点)为同轴圆弧，每个圆弧的圆心都位于中央经线上；中央经线的投影为一直线，其他经线的投影为曲线，并与中央经线对称，且离中央经线越远，其经线的间隔越成比例地递减；极点投影长度为赤道投影长度的一半。
由于等差分纬线多圆锥投影没有给出直接的投影公式，所以用户在使用该投影时往往需要自己来设置一些参考点才能够进行曲线拟合[3]，求解出该投影的正反解变换公式。
MATLAB是一种集数学计算功能、图形化显示功能等多种功能于一体的高级科学计算软件，能够高效快速地解决各种数学相关的科学问题，因此MATLAB被广泛应用于不同领域[4]。本文将在MATLAB平台下进行等差分纬线多圆锥投影的正反解变换，实现该投影下世界地图的矢量和栅格数据输出，以使这些数据可以直接在地理信息系统软件ArcMap中读取和编辑，方便用户直接使用。
1 等差分纬线多圆锥投影变换在MATLAB中的实现
1.1 投影变换使用的地图源数据
(1)世界地图的.shp格式文件(见图1)。
(2)像元行列比为1∶2的栅格世界地图(8192×4096像素)(见图2)。
1.2 投影变换使用的函数文件
以下均为在MATLAB软件平台下进行投影变换所用到的公式。这些公式可以在MATLAB软件中各自存储成一个.m类型的函数文件，以方便其他文件调用。
(1)正解变换公式：(xi，yi)=map_forwar d(long，lati)
式中，long为经度，范围在-180～180之间；lati为纬度，范围在-90-90之间；在经度值为long和纬度值为lati的条件下，xi，yi分别为经等差分纬线多圆锥投影正解变换[1，3，5，6]求得的x值和y值。若long和lati的值超过一定范围时，xi和yi均返回NaN值。
(2)反解变换公式：(long，lati)=map_inverse(xi，yi)
式中各个变量的定义和前面第2节(1)相同，在x、y坐标值分别为xi和yi的条件下long和lati分别为经等差分纬线多圆锥投影反解变换[1，6]求得的经度值和纬度值。若xi和yi的值超过投影区域范围时，long和lati均返回NaN值。
(3)经纬度转像元行列号的公式：(x，y)=image_position(long，lati，image)
式中long和lati的定义和前面第2节(1)相同，image为读入的图像矩阵，在该函数的内部读取该图像的总行数imx和总列数imy。像元行列号x和y分别与long和lati的换算关系如下：
(1)
(2)
式中，[x]代表不超过x的最大正整数。
该式通过指定的经纬度转换成栅格数据所对应的像元的行列号，把该行列号所在的像元值赋值给指定区域。因此该式主要用于栅格数据等差分纬线多圆锥投影的变换。
1.3 投影变换在MATLAB中的实现
(1)矢量数据的投影变换。
以前面所述的.shp矢量世界地图为例，设使用的矢量世界地图文件名为“world _map.shp”，读取、转换、显示和输出中央经线为东经150°的等差分纬线多圆锥投影世界地图代码如下[7]：
S=shaperead('world_map.shp')； %读取shapefile文件。
Slength=length(S)； %读取属性长度。
for cou=1：Slength
xlength=length(S(cou，1).X)；
for count=1：xlength
long=S(cou，1).X(count)；lati=S(cou，1).Y(count)； %读取某一属性的经纬度坐标。
if long
long=long+210； %的经度变换。
else
long=long-150；
end
[xi，yi]=map_forward(long，lati)； %遍历该Shapefile文件上的每一点，通过逐
S(cou，1).X(count)=yi； %点变换，最终生成该投影下的地图数据。
S(cou，1).Y(count)=xi；
end
end
mapshow(S)； %显示投影变换后的地图数据。
shapewrite(S，'world_map_conversion.shp')；%导出地图数据。
在以上代码中，第二行执行后的S变量实际上是一个数组，且数组中的每一个元素都是结构体，上面显示了几何类型、属性名称、X坐标和Y坐标等信息。通过嵌套循环语句逐个访问了每一个点的经纬度坐标，并将变换后的X、Y坐标逐个赋值给S中原来的经纬度坐标。变换后还输出并保存为Shapefile格式文件。变换后的结果如图3所示。
(2)栅格数据的投影变换。
以前面所述的栅格世界地图数据为例，设使用的栅格世界地图文件名为“world_map.jpg”，输出分辨率为5000×4000的中央经线为东经150°的等差分纬线多圆锥投影世界地图。实现这一操作的代码如下：
Mapc=uint8(ones(4000，5000，3))； %建立三维矩阵，图像矩阵类型为uint8。
imagec=imread('world_map.jpg')； %读取栅格世界地图文件。
Mapc=Mapc*255； %使Mapc输出的图片背景为白色。
[p，q，r]=size(Mapc)； %读取Mapc矩阵的尺寸。
for i=1：p
xi=(i-2000)*7000； %第6行和第8行：通过一系列变换，使xi和yi
for j=1：q %的范围能包含投影范围。
yi=(j-2500)*7000；
[long，lati]=map_inverse(xi，yi)； %投影反解变换。
if isnan(long)&&isnan(lati) %若long和lati的值为NaN，不进行颜色赋值，
continue； %并且直接跳入下一轮循环。
end
if long=-180 %此条件语句为输出中央经线为东经150的投影。
long=long+150；
else
if long
long=long-210；
end
end
[row，column]=image_position(long，lati，imagec)； %读取经纬度对应栅格世界
for k=1：r %地图的行列号。
Mapc(i，j，k)=imagec(row，column，k)； %把世界地图的任意一点的
end %像元值赋值给投影地图的
end %对应位置上。
end
image(Mapc)； %显示地图图像。
该代码首先是建立5000×4000的空白图像，然后通过嵌套循环语句逐点扫描该空白图像，把每一个像元的行列号(i，j)经一系列的变换转化成该投影下的x、y坐标，利用反解变换求解对应的经纬度坐标，之后读取经纬度坐标所对应的栅格世界地图图像所对应的行列号(row，column)，把世界地图中行列号为(row，column)的像元值赋值给5000×4000图像的行列号为(i，j)的像元。当循环语句结束后，显示投影变换后的图像并保存。由于保存后的图像已缩小7000倍，这样如果在ArcMap中直接打开该图像，会导致数据显示不正确。但“GeoTiff”解决了这一瓶颈。tfw文件是TIFF文件坐标信息的文本文件，ArcInfo、MicroStation、AutoCAD等软件均支持该格式的坐标信息文件。此文件定义了栅格图像素坐标与实际大地坐标的仿射关系[8，9]。一个.tfw格式的文件构成如下[10]：
第1行：A：地理x坐标中的像元分辨率
第2行：D：y轴旋转系数
第3行：B：x轴旋转系数
第4行：E：地理y坐标中的像元分辨率
第5行：C：左上角第1个像元中心的x坐标值
第6行：F：左上角第1个像元中心的y坐标值
各个参数的图解如图4所示。
由此保存分辨率为5000×4000的投影后的世界地图为TIF格式，用记事本编写对应的tfw格式文件。最后在ArcMap中加载这个世界地图，如图5所示。
2 结语
本文依托MATLAB软件在图形处理方面的优势，综合利用等差分纬线多圆锥投影的正反解变换公式，还利用了像素长宽比为2∶1的世界地图原始数据中经纬度与行列号的关系，实现了矢量数据和栅格数据的等差分纬线多圆锥投影变换。此外本文以桌面GIS软件ArcMap为例，使用了.tfw格式文件让输出后栅格数据等差分纬线多圆锥投影世界地图准确显示。本文中的投影变换呈现出如下特点。
(1)输出后的矢量和栅格地图数据可以供用户在ArcMap、MapInfo等GIS软件中被直接加载，方便用户进行编辑。
(2)输出后的矢量和栅格地图数据具有准确性、可靠性。与文献[3]相比，整个变换过程无需人工干预，节省了人工计算配准点的繁重劳动，减少了栅格地图数据地理配准中可能出现的误差。
参考文献
[1] 董曼，李胜乐.世界地图等差分纬线多圆锥投影的正反解变换[J].大地测量与地球动力学，2008，28(2)：95-99.
[2] 祝国瑞.地图学[M].武汉：武汉大学出版社，2004.
[3] 叶远智.栅格数据的等差分纬线多圆锥投影转换[J].测绘通报，2012(9)：68-70.
[4] 刘正君.MATLAB科学计算宝典[M].北京：电子工业出版社，2012.
[5] 胡毓钜.地图投影[M].北京：测绘出版社，1981.
[6] 杨启和.地图投影变换原理和方法[M].北京：解放军出版社，1989.
[7] 苏金明，王永利.MATLAB图形图像[M].北京：电子工业出版社，2005.
[8] 余建军，蔡以雷，洪景峰.一种利用数字线划图快速生产数字栅格地图的方法[J].测绘通报，2009(7)：51-54.
[9] 王宪民.ArcGIS在栅格数据处理中的应用初探[J].测绘技术装备，2005，1(7)： 32-34.
[10] 吴铭杰.基于AutoCAD的TIFF world(TFW)图像的应用研究[A].史照良.第十四届华东六省一市测绘学会学术交流会论文集[C].南京：《现代测绘》编辑部，2012：157-158.
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
