地图投影设置

在上边我们已经用到了一些投影参数。比如进行投影的反转,所使的是-I

中国等积投影是:

proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25

其中有许多参数前边都加了前缀+, 这后边的参数是对地图投影的真正的设置,proj命令也将按这个规定来进行转换。 这种的参数的形式是 +param=valueparam(参数)来给定一个 value(值), 这个值可以是一个以度分秒的格式也可以是实数,整数格式,甚至可以是一个 ASCII 字符串。 需要注意的是一种拼错了的参数名不会被理会,如果一个参数输入了两次, 则第一次输入的将会被使用。所以在执行前要将参数检查好。 另一个很有用的特点是, proj会自动决定中央经线并且追加一个+lon_0参数到你的定义中。

选择投影

选择投影使用的参数是+proj=name。 proj中集成了不少的投影而且还在不断的增加。 可惜的是我们中国的投影并没增加多少。 在这里提供我们自己增加投影的方法。我们可以用下面命令:

$proj -l #得到内置的投影类型
aea : Albers Equal Area #可以看到 Albers等积投影的在 proj中的值是 aea
$proj -lP #得到更详细的投影信息,投影包含的参数。
aea : Albers Equal Area   #更详细的信息
Conic Sph&Ell 
lat_1= lat_2= 

这样我们就可以定义你想要的投影了。

选择椭球体

虽然投影已经选定了,可是还要确定椭球体。有两种方法:一种是利用内置的一些椭球体。

$proj -le #显示出内置的椭球体。

以下是我国经常用的椭球体:

clrk66 a=6378206.4      b=6356583.8      Clarke 1866 
krass a=6378245.0      rf=298.3         Krassovsky, 1942 
WGS84 a=6378137.0      rf=298.257223563 WGS 84

另一种就是自定义一个你自己的地球。

自定义椭球体参数

上面的三个椭球体上就包含了定义一个椭球体的参数。 下面是这几个参数的含义:必须有的:

+a=A   椭球的赤道半径(半长轴)

下边的这些参数只需要一个就可以:

  • +b=B 椭球的极半径(半短轴)
  • +f=F 椭球的扁率 $F=\frac{(A-B)}{A}$
  • +rf=RF 椭球的反扁率 $RF=\frac{1}{F}$
  • +e=E 偏心率
  • +es=ES 偏心率的平方 $\frac{E}{2}$

例如指定一个 Clark1866 椭球来作为中国等积投影的参数,可以这样设置:

$proj +proj=aea +ellps=clrk66 +lon_0=105 +lat_1=25 
105 36                  #可以和上边的 krass椭球体比较
0.00 3850517.66 
104 36 
-88731.89 3850957.52 
106 24 
102046.72 2507997.23

同样也可以自定义:

$ proj +proj=aea +a=6378206.4 +es=.006768658 +lon_0=105 +lat_1=25

或:

$ proj +proj=aea +a=6378206.4 b=6356583.8 +lon_0=105 +lat_1=25

得到结果都是一样的。如果只指定了一个+a则会得到一个正规的球体。例如:

proj +a=1

你会定义一个半径 R = 1的单位球体。

其它的参数

我们知道 UTM投影,为确保每个带中的点全为正, 在北半球将坐标从每带的中央经线西移 500公里。而南半球为了保证 Y轴为正,不得不还向南移。 这样两个参数在 proj中可以用两个参数+x_0+y_0来确定。例如在UTM中,

+x_0=5000000 +y_0=0

还有一个是+lat_0

第四个参数lat_0 = 0用于指定几个投影的中心平行和相关的y轴原点。

还有最后一个+lon_0:指中央经线,在 UTM和 Albers中要设置。 上边的四个参数是可选的,如果没有将会默认为 0,除+lon_0,它会计算。 这样我们就完成了一个投影的完整定义了, 我认为这是一种很简单很优美的投影定义方式。通常把这种 proj 定义出的投影形式称为 PROJ.4 格式,例如中国的等积投影:

+proj=aea +ellps=krass +lon_0=105 +lat_1=25

还有中国的 Lambert投影,写法和上边差不多,在 proj4里的投影名称为 +proj=lcc

一些实例

定义 UTM 投影的方法:

$proj +proj=utm +zone=48  
105 36   #在离兰州最近的一条中央经线
500000.00 3983948.45  #可以看出西移 500000米

定义高斯克吕格(它是横轴墨卡托的一种,也是横轴黑卡托的默认值)投影的方法:

$proj +proj=tmerc +lon_0=105  
#这上边没带,可以简单的算出每个带的中央经线的度数:(带数*6)-3就是中央经线

自定义你的 proj 参数

我们希望将我们常用的地图投影的定义放到文件中,我们只需在需要时将其导入即可。 PROJ.4 也提供了这个功能。我们以设置+init参数来进行。 首先先将地图投影保存到文件中。我们保存地图投影的文件格式如下: proj.dat

在安装 PROJ.4 的位置,有一个默认的文件夹,来指示不满足要求的初始化文件名称, 以及投影的默认文件 proj_def.dat。 环境变量 PROJ_LIB 可以用来给出另一个新的文件夹。

cs2cs

cs2cs程序允许任何坐标系统之间的转换,支持不同基准面的转换。 cs2cs根据一些输入的点进行原坐标系统和目标坐标系统之间的转换。 既能够进行地理坐标与投影坐标的转换,又能实现不同基准面之间的转换。

下面是命令的用法:

cs2cs [ -eEfIlrstvwW [ args ] ] [ +opts[=arg] ]
   [ +to [+opts[=arg]] ] file[s]

下面是cs2cs的参数说明,且控制参数可能以任意形式出现:

| ll 参数 & 说明
| -I & 指定逆翻译的方法,从坐标系统到原坐标系统
| -ta & A specifies a character employed
| -e string &
字符串是任意字符串,如果在数据转换过程中发现错误的话。默认值是t.
| -E & 复制输入坐标
| -l & 投影标识符
| -r & 反转输入顺序
| -s & 反转输出顺序
| -f format & 输出格式
| -[w\|W]n & 保留小数位数
| -v & 对制图控制参数进行测试,不与——T选项一起使用

运行行参数与制图参数相互联系,用法随着投影的改变而改变。

cs2cs程序需要两个投影系统的定义。第一个坐标系统的定义,是根据没有在 +to argument后出现的投影参数, 所有在 +to argument后的投影参数都被看作是第二个投影系统的定义。 如果没有第二个系统的定义,那么就会假定与原坐标系统的基准和椭球相同。 原系统和新系统都是地理坐标系统,都可以进行投影,但可能它们之间有不同的基准面。

额外的投影控制参数可能包含在两个附加的控制文件中: 第一个文件以+init=file:id为参考,第二个文件是在投影名称建立之后处理。 参数PROJ_LIB创建了文件的默认目录。这也常用来支持基准面转换文件。

一个或者多个文件需要转换,需要按照从左到右的顺序进行转换。 A会确定处理的标准和输入的位置。如果没有文件指定,那么输入将从stdin开始。 对于输入数据,两个数据值必须是前两个空格分隔的字段。

输入的地理数据必须是DMS格式, 输入的笛卡尔数据单元必须与椭球的主轴和所求半径单元一致。 输出的地理坐标为DMS格式,且输出值精确到0.001。零值的分秒数据被删除。

下面是一个实例:

cs2cs +proj=latlong +datum=NAD83 +to +proj=utm +zone=10  +datum=NAD27 -r <<EOF
45d15'33.1"   111.5W
45d15.551666667N   -111d30
+45.25919444444    111d30'000w
EOF

将会转变输入的NAD83地理坐标标成NAD27,投影为UTM10带。 这个例子的地理值是等值的,是DMS输入众多形式的一个例子。

1402285.99      5076292.42 0.000

geod

geod(直接)和invgeod(间接)执行测地线(大圆)计算, 用于确定给定初始点的纬度、经度、方位角和距离(直接)或前后方位角的终点的纬度、经度和后方位和初始纬度和终点纬度之间的距离以及距离(间接)。

geod 有两个程序。

  • geod - 直接测地计算
  • invgeod - 间接测地计算

下面是geod和invgeod用法:

geod  +ellps=<ellipse> [ -afFIlptwW [ args ] ] [+args ] file[s]
invgeod +ellps=<ellipse> [ -afFIlptwW [ args ] ] [ +args ] file[s]

下面是程序运行的参数说明,且运行行控制参数能够以任何顺序出现:

ll 参数 & 说明
-I & 指定逆短程计算
-a & 初始点和终点的经度和纬度
-ta & 表示不处理的话控制线就会通过
-le & 给出所有椭球的列表
-lu & 给出所有单位的列表
-[f|F] & 格式
-[w|W]n & 小数点保留位数
-p & 产生方位值,0-360之间

运行行参数和大地参数一起用于定义椭球或者是球。看proj文档并获取这些参数的完整列表。选项从左到右运行,如果第一个出现值被认为是期望的值,则重新进入命令是忽略的。

一个或者多个文件需要转换,按照从左到右的顺序。A会确定处理的标准输入的位置。如果没有文件指定,那么输入将从stdin开始。

对于直接输入数据,必须按照纬度,经度,方位和距离的顺序,输出依次是纬度,经度和总站点方位。初始和总站点的纬度经度是逆模式的输入,初始点和总站点的前向,逆向方位和点与点之间的距离是输出。

输入的地理坐标和方位数据必须是DMS格式,输入的距离数据必须与椭球的主轴和球的半径单位一致。输出地理坐标是DMS格式,输出的距离数据与椭球和球的半径单位相同。

地球椭球的选择和程序与proj的选择是一样的。

GEOD可以用来确定沿两条地理线之间的中心点,也可以确定沿弧线的中心点。 在这两种情况下,初始点必须被指定,用 +lat_1=latlon_1=lon为参数。 总站点为+lat_2=lat+lon_2=lon,到初始点的距离和方位这样指定:+S=distance and +A=azimuth

如果点是按照测地线确定,那么不是用+n_S=integer来确定中点的个数就是用+del_S=distance来指定点之间的距离。

沿弧线确定点的话,+del_A=angle+n_A=integer都得确定,他们决定了各自的角增量和点的数量。

以下脚本确定了从美国马萨诸塞州波士顿到波特兰的美国法定里程的测地方位和距离:

geod +ellps=clrk66 <<EOF -I +units=us-mi
42d15'N 71d07'W 45d31'N 123d41'W
EOF

得到了以下结果:

-66d31'50.141"   75d39'13.083"   2587.504

其中前两个值是从波士顿到波特兰的方位角,而从波特兰到波士顿的后方位角跟随距离的。

前向测地线使用的例子是使用波士顿的位置,并以方位角和距离来确定波特兰的位置:

geod +ellps=clrk66 <<EOF +units=us-mi
42d15'N 71d07'W -66d31'50.141" 2587.504
EOF

结果:

45d31'0.003"N    123d40'59.985"W
75d39'13.094"

注意:距离值的精度并不会影响波特兰位置的精度。

proj命令的使用

PROJ.4 里集成了许多制作地图用的投影参数,并且实现了一些命令行程序,从而进行地图投影、投影转换,以及其他投影相关的操作。

  • proj 程序仅限于同一基准的地理坐标与投影坐标之间的转换。
  • invproj 逆向制图投影转换,即 proj 的即变换。
  • cs2cs 与上面两个程序类似,而且允许任意地图投影坐标系统之间的转换,包括不同基准面之间的转换。
  • nad2nad 程序提供NAD27与NAD83之间的转换(同样的功能在 cs2cs 中已经实现了,但是这个更加方便)。
  • geod 程序提供了大圈计算的功能。

proj 是投影处理的主要的命令,我们重点来看一下此命令的用法。 需要注意的是,proj 在进行投影转换时,是使用交互方式运行的。 当使用命令行设定的正确的投影空间参考环境后。交互式输入经纬度坐标, 会产生在此空间参考下的投影后的空间坐标。

先来了解一下 proj 命令的各个参数。

投影转换

我国常用的地图投影包括Albers,Lambert,Gauss-Kruger,UTM投影等。

等积投影由于没有面积变形, 所以在土地调查,植被盖度分类等涉及到要保持面积不能变形的情况下, 中国的全国性地图大多采用等积投影。国际上称等积投影为 Albers 投影,这是一种圆锥等积投影。中国所使用的 Albers 的参数是双标准纬线, 25N , 47N , 中央经线为 105E , 椭球体为Krassovsky。用 PROJ.4 表示为:

+proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 

下边将使用中国的 Albers投影,简称为 Albers_China,来作个简单的投影转换。 命令的作用是将经纬度坐标转换为定义好的坐标系统中的坐标。

$proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 
105 36
0.00 3847866.97 
104d30' 36d30'
-43977.163904491.79
104.5N 36.5E
-43977.163904491.79

上面命令中的第1,2,4,6都是输入的命令,3,5,7是程序返回的结果。 第一行说明要进行投影转换(经纬度至地理坐标),并给出了相关的参数。 然后下面每一行输入两个数值,数值之间使用空格分隔。 数据输入可以用度分秒表示,也可以用十进制的度表示,并且可以指定南北纬与东西经。

同样也可以进行反转,即将Albers转为经纬度,只要在命令中加入参数 -I

$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 -I
102064.08 2503934.26 
106dE  24dN 

经纬度的反转输入与输出

在这里转换的过程中始终是按经度、纬度 (x, y) 的顺序放进的。

如果输入时想掉转可在命令中加 -r ,如果是输出想掉转,可以在命令中加 -s

$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 -r -s 
36 105 
3847866.97  0.00 
33 104 
3509623.92  -91933.97 

命令行中的EOF

Proj.4 是一个典型的Unix程序,除了最基本的用法,还可以使用Unix的管道工具, 从而来简化(加强)输入与输出。

例如,可以进行批量转换:

$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 <<EOF 
> 105 36 
> 104 36 
> 106 24 
> EOF 
0.00 3847866.97 
-88522.43 3848312.80 
102064.08 2503934.26 

使用文件来进行处理

Proj.4同样也可以使用文件来进行批量转换:

对于文本文件 lat_lon.test ,其内容如下:

105dE 36dN 
104dE 36dN 
106dE 24dN 

使用下面的命令进行转换:

$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 ~/lat_lon.test > alberst.test 

生成的: alberst.test ,内容如下:

0.00 3847866.97 
-88522.43 3848312.80 
102064.08 2503934.26 

同时你也可以在文件中加注释和对坐标点进行说明,在转换后仍可以保留:

lat_lon.test 
#it's just a test for convert file format 
105dE 36dN not Lanzhou 
104dE 36dN Lanzhou 
106dE 24dN Unknow place 

命令:

$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 
+lat_2=47 ~/lat_lon.test >albers.test 

生成的albers.test,内容如下:

#it's just a test for convert file format 
0.00 3847866.97 not Lanzhou 
-88522.43 3848312.80 Lanzhou 
102064.08 2503934.26 Unknow place 

在命令上边的~/lat_lon.test是输入的文件在下指的是当前目录, Windows下没试过,不过可以用绝对路径。 >是重定向,输出文件。

地图单位

PROJ.4 的默认单位为米(meter), 我们设置参数 +units 来控制输入的坐标单位。 我们可以将它的输入或输出的数据的单位改为其它:

$ proj +proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47 +units=km -I <<EOF
> 0 3847.86697 
> -88.52243 3848.31280 
> 102.06408 2503.93426 
> EOF
105dE  36dN 
104dE  36dN 
106dE  24dN