MATLAB地图工具箱学习总结（三）地图工具箱的基本知识

MATLAB地图工具箱学习总结系列：

（一）从地图投影说起

（二）大圆和恒向线

（三）地图工具箱的基本知识

（四）自定义投影

今天想要介绍的是一些比较基础的函数。了解了这些函数，地图投影的基本概念才能真正明白。而要想继续研究MATLAB中有关地图投影的函数，尤其是未来我要提到的投影文件源代码，知晓这些函数的功能必不可少。本篇文章将会罗列三个案例，并在后面一一进行讲解。

1                    作业案例：地图投影作业1

这次的案例从作业1开始。作业1是要求计算出地球椭球体的一些基本参数，包括子午圈曲率半径、卯酉圈曲率半径、平均曲率半径和纬圈半径等。当初我交上的作业完全是数学公式的堆砌，不过其实MATLAB中有相关的函数。

在这里用到了三个函数，包括referenceEllipsoid, rcurve和rsphere，其中angle是输入的角度。利用这两个函数轻松就将以上四种椭球体的参数获得了。

让我们来分析一下代码。

首先需要介绍的是referenceEllipsoid（参考椭球体）函数。和这个函数相似的还有referenceSphere（参考球体）和oblateSpheroid函数，都是关于参考系的设置。函数括号中的参数也很明显，分别是设置为WGS84坐标系，单位是km.在进行地图投影的计算前，设置好所需的参考椭球是必须的。MATLAB中包含了多种参考椭球，可以直接使用，只需要查看函数帮助就可以了。

接下来要讲得是rcurve函数。这个函数可以计算各种曲率半径。其基本用法为：rcurve(propertyName,ellipsoid,lat).在上面案例中提到的三个参数分别是：transverse，用来计算卯酉圈曲率半径；meridian，用来计算子午圈曲率半径；parallel，用来计算纬线圈半径。只要接下来输入所用参考椭球和纬度，即可计算出相应的参数。

最后要讲一下rsphere函数。这个函数用来计算地球球体的半径。它有多种计算方式，其中triaxial参数可以计算出平均曲率半径。公式为sqrt(a*b),即几何平均值。此外还包括authalic、euler等参数，具体形式请查看帮助。

2                    作业案例2：地图投影作业2

作业2的问题是，一个人从哪里出发，先向东走，再向北走（过了北极后方向仍然不变），最终回到原点，且向东向北走的距离相同。其核心问题是求经线和纬线长度。还记得当时把一整长串积分公式好不容易录入MATLAB，MATLAB向我报告积分太多算不了。后来改用了其泰勒展开式的一部分带入进行计算，好在成功了。不过在本文里我并没有更好的解法。但我想借此机会介绍几个函数。分别是departure,meridianarc,meridianfwd,distance,reckon。

departure函数的作用在于计算纬线的长度。其基本用法是departure(lon1,lon2,lat,ellipsoid)，首先输入两条经线，再输入纬度，带入椭球体参数即可计算出相应纬线的长度。

meridianarc、meridianfwd函数刚好与departure函数相反，是用来求经线长度的。meridianarc的基本用法是s=meridianarc(phi1,phi2,ellipsoid)，即只要输入纬度即可算出它们间相隔的距离s。meridianfwd则和meridianarc互为正反算。其基本用法是phi2=meridianfwd(phi1,s,ellipsoid)，即输入一点经度，再输入相隔的距离和椭球体，即可算出相对应的经度。需要注意到的是，这两个函数都需要输入弧度制的经度，所以一定要先转换好再输入。

以上三个函数的效果如下：

distance和reckon两个函数都和距离有关。distance的基本用法是[arclen,az]=distance(lat1,lon1,lat2,lon2),很明显，只要输入两个点的经纬度，即可获得两点间的距离。这个距离默认为大圆距离，当然也可以设置rh，求恒向线的距离。获得的arclen即为距离，此外还可以获得az，即两点之间的方位角。reckon的基本用法是[latout,lonout]=reckon(lat,lon,arclen,az)，这和distance刚好相反。输入一个点的经纬度，距离和方位角，即可求出距离这个点相应距离和相应方位角的点。

3                    作业案例3：地图投影作业7

在作业7中，老师要求我们每个人根据标准纬线，使用圆锥等角投影显示各自家乡的省份。我想在这里提一提MATLAB中地图投影最基本的组织原理。

在工具箱中，有一个mstruct的句柄，这便是地图投影组织结构。建立axesm后，使用命令：gcm就可以获得当前的投影信息。我们以墨卡托投影为例来看一下里面包含了哪些结构。

当然不仅仅这些，下面还有很多属性没有列出。这里面就显示了axesm的各种属性。而这些属性我在本系列的第一篇中就提到了，虽然当时的我还没有理解其中的原理。

对于属性的设置，除了用axesm可以进行设置外，还可以分别使用getm和setm对属性进行读取和设置。getm就是读取投影属性，使用方法为：getm(gcm,propertyName),即输入你想要了解的属性即可。而setm则需要在propertyName后加具体的值，那么地图投影就会跟着改变。

其中有一个mapparallels是解决这道题的关键，这个属性正是指的地图投影的标准纬线。nparallels可以设定标准纬线的个数，当nparallels为1时，可以设定一个mapparallels，说明是切投影，两个值的时候则是割投影。要注意的是，并不是所有的投影类型都支持设置该属性，因此在使用前要确认这个投影是否能设置标准纬线。

不知不觉，当我写到这个时候，地图投影的课程已经结课了。最后的作业是自定义投影，也是花了我不少时间，不过还是做出来了。那么，下一篇也就将是本系列的最后一篇了，希望能够给使用MATLAB地图工具箱的人带来更多的帮助。

天靖居士

2016.5.16

8.17更新说明：具体代码请参考：https://git.oschina.net/kkyyhh96/MapProjectInMatlab

有关地图工具箱其他文章，请参看：

MATLAB地图工具箱学习总结系列：

（一）从地图投影说起

（二）大圆和恒向线

（三）地图工具箱的基本知识

（四）自定义投影