>>> from env_helper import info; info()
页面更新时间: 2023-04-15 20:04:30
运行环境:
    Linux发行版本: Debian GNU/Linux 12 (bookworm)
    操作系统内核: Linux-6.1.0-7-amd64-x86_64-with-glibc2.36
    Python版本: 3.11.2

6.1. PROJ.4简介

PROJ.4是开源GIS最著名的地图投影库,它专注于地图投影的表达,以及转换,许多GIS开源软件的投影都直接使用Proj.4的库文件。GDAL中的投影转换函数(类CoordinateTransformation中的函数)也是动态调用该库函数的。

Proj.4的功能主要有经纬度坐标与地理坐标的转换,坐标系的转换,包括基准变换等。地图投影的表达方式有多种,由于采用一种非常简单明了的投影表达--PROJ.4比其它的投影定义简单,很容易就能看到各种地理坐标系和地图投影的参数,同时它强大的投影转换功能,也是非常吸引人的。

PROJ.4 在 Window 的命令下有可运行的 EXE文件,其实它更主要的是一个类库!在 Linux下除了可以直接运行外,还可以作为库来进行更高功能的开发。可以用来开发一些批处理功能。

6.1.1. 安装

Window下安装

从 PROJ.4 的网站下载,解压缩,并安装文件(bin压缩包)。 按照Readme说明用DOS命令将文件夹添加到系统环境变量, 之后可用示例数据测试经纬度与平面坐标的转换。(DOS命令不会的自行解决。)

PROJ.4的官方网站为 http://trac.osgeo.org/proj/

用DOS命令添加系统环境变量,当此次命令结束,系统环境变量即恢复原状,下次使用仍需要添加。 若在“我的电脑”属性里添加系统环境变量,则一次设置可多次使用。若进行开发的话,可以使用pyproj。

Debian/Ubuntu Linux下安装

在Debian Stretch中,可以在终端输入:

# apt-get install proj-bin python3-pyproj

其中proj-bin为 PROJ.4 的命令行工具,实现了在Linux下的Shell环境中进行地图投影处理的功能; python3-pyproj则是在 Python 3中的实现。

pyproj在 Windows 和 Linux下都很好安装。 pyproj 原来托管在 Google Code,后来迁移到 GitHub https://github.com/jswhit/pyproj

Pyporj是 Python下的 proj。 可以很方便的对点进行地图投影转换。 同时还可以在它的基础上开发出更高级的应用。 Pyproj包里包括两个类, Proj 类和 Geod 类。

安装完成后,使用下面语句导入Proj类

pip安装

pip环境下的安装,使用如下:

pip install pyproj

运行下面代码,查看安装是否成功:

>>> from pyproj import Proj

6.1.2. 在PROJ.4中了解基准面与椭球体

使用proj命令查看椭球体参数

proj 命令是对地球经纬度进行投影的,即将经纬度坐标转换为地理坐标。 当然也可以将地理坐标转换为经纬度坐标。即在终端下输入:

>>> !proj
Rel. 9.1.1, December 1st, 2022
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]

会显示出 proj 程序的用法。包括参数设置,可选项,和输入文件。

6.2. 显示参数

我们可以使用下边的命令来显示 proj 里内置的有关地图投影的参数。 显示投影类型:

>>> !proj -l | wc -l
178
>>> !proj -l
>>>
adams_hemi : Adams Hemisphere in a Square
adams_ws1 : Adams World in a Square I
adams_ws2 : Adams World in a Square II
aea : Albers Equal Area
aeqd : Azimuthal Equidistant
affine : Affine transformation
airy : Airy
aitoff : Aitoff
alsk : Modified Stereographic of Alaska
apian : Apian Globular I
august : August Epicycloidal
axisswap : Axis ordering
bacon : Bacon Globular
bertin1953 : Bertin 1953
bipc : Bipolar conic of western hemisphere
boggs : Boggs Eumorphic
bonne : Bonne (Werner lat_1=90)
calcofi : Cal Coop Ocean Fish Invest Lines/Stations
cart : Geodetic/cartesian conversions
cass : Cassini
cc : Central Cylindrical
ccon : Central Conic
cea : Equal Area Cylindrical
chamb : Chamberlin Trimetric
collg : Collignon
col_urban : Colombia Urban
comill : Compact Miller
crast : Craster Parabolic (Putnins P4)
defmodel : Deformation model
deformation : Kinematic grid shift
denoy : Denoyer Semi-Elliptical
eck1 : Eckert I
eck2 : Eckert II
eck3 : Eckert III
eck4 : Eckert IV
eck5 : Eckert V
eck6 : Eckert VI
eqearth : Equal Earth
eqc : Equidistant Cylindrical (Plate Carree)
eqdc : Equidistant Conic
euler : Euler
etmerc : Extended Transverse Mercator
fahey : Fahey
fouc : Foucaut
fouc_s : Foucaut Sinusoidal
gall : Gall (Gall Stereographic)
geoc : Geocentric Latitude
geogoffset : Geographic Offset
geos : Geostationary Satellite View
gins8 : Ginsburg VIII (TsNIIGAiK)
gn_sinu : General Sinusoidal Series
gnom : Gnomonic
goode : Goode Homolosine
gs48 : Modified Stereographic of 48 U.S.
gs50 : Modified Stereographic of 50 U.S.
guyou : Guyou
hammer : Hammer & Eckert-Greifendorff
hatano : Hatano Asymmetrical Equal Area
healpix : HEALPix
rhealpix : rHEALPix
helmert : 3(6)-, 4(8)- and 7(14)-parameter Helmert shift
hgridshift : Horizontal grid shift
horner : Horner polynomial evaluation
igh : Interrupted Goode Homolosine
igh_o : Interrupted Goode Homolosine Oceanic View
imoll : Interrupted Mollweide
imoll_o : Interrupted Mollweide Oceanic View
imw_p : International Map of the World Polyconic
isea : Icosahedral Snyder Equal Area
kav5 : Kavrayskiy V
kav7 : Kavrayskiy VII
krovak : Krovak
labrd : Laborde
laea : Lambert Azimuthal Equal Area
lagrng : Lagrange
larr : Larrivee
lask : Laskowski
lonlat : Lat/long (Geodetic)
latlon : Lat/long (Geodetic alias)
lcc : Lambert Conformal Conic
lcca : Lambert Conformal Conic Alternative
leac : Lambert Equal Area Conic
lee_os : Lee Oblated Stereographic
loxim : Loximuthal
lsat : Space oblique for LANDSAT
mbt_s : McBryde-Thomas Flat-Polar Sine (No. 1)
mbt_fps : McBryde-Thomas Flat-Pole Sine (No. 2)
mbtfpp : McBride-Thomas Flat-Polar Parabolic
mbtfpq : McBryde-Thomas Flat-Polar Quartic
mbtfps : McBryde-Thomas Flat-Polar Sinusoidal
merc : Mercator
mil_os : Miller Oblated Stereographic
mill : Miller Cylindrical
misrsom : Space oblique for MISR
moll : Mollweide
molobadekas : Molodensky-Badekas transformation
molodensky : Molodensky transform
murd1 : Murdoch I
murd2 : Murdoch II
murd3 : Murdoch III
natearth : Natural Earth
natearth2 : Natural Earth 2
nell : Nell
nell_h : Nell-Hammer
nicol : Nicolosi Globular
nsper : Near-sided perspective
nzmg : New Zealand Map Grid
noop : No operation
ob_tran : General Oblique Transformation
ocea : Oblique Cylindrical Equal Area
oea : Oblated Equal Area
omerc : Oblique Mercator
ortel : Ortelius Oval
ortho : Orthographic
pconic : Perspective Conic
patterson : Patterson Cylindrical
peirce_q : Peirce Quincuncial
pipeline : Transformation pipeline manager
poly : Polyconic (American)
pop : Retrieve coordinate value from pipeline stack
push : Save coordinate value on pipeline stack
putp1 : Putnins P1
putp2 : Putnins P2
putp3 : Putnins P3
putp3p : Putnins P3'
putp4p : Putnins P4'
putp5 : Putnins P5
putp5p : Putnins P5'
putp6 : Putnins P6
putp6p : Putnins P6'
qua_aut : Quartic Authalic
qsc : Quadrilateralized Spherical Cube
robin : Robinson
rouss : Roussilhe Stereographic
rpoly : Rectangular Polyconic
s2 : S2
sch : Spherical Cross-track Height
set : Set coordinate value
sinu : Sinusoidal (Sanson-Flamsteed)
somerc : Swiss. Obl. Mercator
stere : Stereographic
sterea : Oblique Stereographic Alternative
gstmerc : Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)
tcc : Transverse Central Cylindrical
tcea : Transverse Cylindrical Equal Area
times : Times
tinshift : Triangulation based transformation
tissot : Tissot
tmerc : Transverse Mercator
tobmerc : Tobler-Mercator
topocentric : Geocentric/Topocentric conversion
tpeqd : Two Point Equidistant
tpers : Tilted perspective
unitconvert : Unit conversion
ups : Universal Polar Stereographic
urm5 : Urmaev V
urmfps : Urmaev Flat-Polar Sinusoidal
utm : Universal Transverse Mercator (UTM)
vandg : van der Grinten (I)
vandg2 : van der Grinten II
vandg3 : van der Grinten III
vandg4 : van der Grinten IV
vertoffset : Vertical Offset and Slope
vitk1 : Vitkovsky I
vgridshift : Vertical grid shift
wag1 : Wagner I (Kavrayskiy VI)
wag2 : Wagner II
wag3 : Wagner III
wag4 : Wagner IV
wag5 : Wagner V
wag6 : Wagner VI
wag7 : Wagner VII
webmerc : Web Mercator / Pseudo Mercator
weren : Werenskiold I
wink1 : Winkel I
wink2 : Winkel II
wintri : Winkel Tripel
xyzgridshift : Geocentric grid shift

在 Debian Wheezy中,一共有投影类型126个,Debian Jessie中为132 个,在Debian Stretch中则为140个。 为缩减篇幅,只列出部分投影。其中第一个就是中国常用的阿尔波斯投影。

PROJ.4支持许多长度单位,可以通过参数 -lu,看到支持的单位:

>>> !proj -lu | wc -l
21
>>> !proj -lu
    mm 0.001                millimetre
    cm 0.01                 centimetre
     m 1                    metre
    ft 0.3048               foot
 us-ft 0.304800609601219    US survey foot
  fath 1.8288               fathom
   kmi 1852                 nautical mile
 us-ch 20.1168402336805     US survey chain
 us-mi 1609.34721869444     US survey mile
    km 1000                 kilometre
ind-ft 0.30479841           Indian foot (1937)
ind-yd 0.91439523           Indian yard (1937)
    mi 1609.344             Statute mile
    yd 0.9144               yard
    ch 20.1168              chain
  link 0.201168             link
    dm 0.01                 decimeter
    in 0.0254               inch
ind-ch 20.11669506          Indian chain
 us-in 0.025400050800101    US survey inch
 us-yd 0.914401828803658    US survey yard

同样的,还有参数 -le ,显示支持的椭球体(ellipsoid)信息, 以及各个椭球体向WGS 84椭球体的转换参数。

>>> !proj -le | wc -l
46
>>> !proj -le
    MERIT a=6378137.0      rf=298.257       MERIT 1983
    SGS85 a=6378136.0      rf=298.257       Soviet Geodetic System 85
    GRS80 a=6378137.0      rf=298.257222101 GRS 1980(IUGG, 1980)
    IAU76 a=6378140.0      rf=298.257       IAU 1976
     airy a=6377563.396    rf=299.3249646   Airy 1830
   APL4.9 a=6378137.0      rf=298.25        Appl. Physics. 1965
    NWL9D a=6378145.0      rf=298.25        Naval Weapons Lab., 1965
 mod_airy a=6377340.189    b=6356034.446    Modified Airy
   andrae a=6377104.43     rf=300.0         Andrae 1876 (Den., Iclnd.)
   danish a=6377019.2563   rf=300.0         Andrae 1876 (Denmark, Iceland)
  aust_SA a=6378160.0      rf=298.25        Australian Natl & S. Amer. 1969
    GRS67 a=6378160.0      rf=298.2471674270 GRS 67(IUGG 1967)
  GSK2011 a=6378136.5      rf=298.2564151   GSK-2011
   bessel a=6377397.155    rf=299.1528128   Bessel 1841
 bess_nam a=6377483.865    rf=299.1528128   Bessel 1841 (Namibia)
   clrk66 a=6378206.4      b=6356583.8      Clarke 1866
   clrk80 a=6378249.145    rf=293.4663      Clarke 1880 mod.
clrk80ign a=6378249.2      rf=293.4660212936269 Clarke 1880 (IGN).
      CPM a=6375738.7      rf=334.29        Comm. des Poids et Mesures 1799
   delmbr a=6376428.       rf=311.5         Delambre 1810 (Belgium)
  engelis a=6378136.05     rf=298.2566      Engelis 1985
  evrst30 a=6377276.345    rf=300.8017      Everest 1830
  evrst48 a=6377304.063    rf=300.8017      Everest 1948
  evrst56 a=6377301.243    rf=300.8017      Everest 1956
  evrst69 a=6377295.664    rf=300.8017      Everest 1969
  evrstSS a=6377298.556    rf=300.8017      Everest (Sabah & Sarawak)
  fschr60 a=6378166.       rf=298.3         Fischer (Mercury Datum) 1960
 fschr60m a=6378155.       rf=298.3         Modified Fischer 1960
  fschr68 a=6378150.       rf=298.3         Fischer 1968
  helmert a=6378200.       rf=298.3         Helmert 1906
    hough a=6378270.0      rf=297.          Hough
     intl a=6378388.0      rf=297.          International 1924 (Hayford 1909, 1910)
    krass a=6378245.0      rf=298.3         Krassovsky, 1942
    kaula a=6378163.       rf=298.24        Kaula 1961
    lerch a=6378139.       rf=298.257       Lerch 1979
    mprts a=6397300.       rf=191.          Maupertius 1738
 new_intl a=6378157.5      b=6356772.2      New International 1967
  plessis a=6376523.       b=6355863.       Plessis 1817 (France)
     PZ90 a=6378136.0      rf=298.25784     PZ-90
   SEasia a=6378155.0      b=6356773.3205   Southeast Asia
  walbeck a=6376896.0      b=6355834.8467   Walbeck
    WGS60 a=6378165.0      rf=298.3         WGS 60
    WGS66 a=6378145.0      rf=298.25        WGS 66
    WGS72 a=6378135.0      rf=298.26        WGS 72
    WGS84 a=6378137.0      rf=298.257223563 WGS 84
   sphere a=6370997.0      b=6370997.0      Normal Sphere (r=6370997)

注意:最后一个,是一个球,而不是椭球。

参数-ld, 显示Proj4支持的基准面(Datum)信息。

>>> !proj -ld
Rel. 9.1.1, December 1st, 2022
<proj>:
invalid list option: ld
program abnormally terminated

可以看到,WGS84是目前最常用的椭球体,且其他椭球体的定义,都是相对WGS84的参数而定义的(towgs84)。

6.2.1. 大地水准面与椭球体

地球是一个表面很复杂的球体, 人们以假想的平均静止的海水面形成的“大地体”为参照, 推求出近似的椭球体,并通过了理论和实践证明, 该椭球体以地球短轴为轴的椭圆而旋转的椭球面, 而且这个椭球面也可以用数学公式来进行表达。

为此,人们从数学角度出发,选择了一个形状和大小与大地球体均极为相似的旋转椭球面来描述地球表层,这也被称为地球椭球面。

大地水准面

众所周知,地球是一个近似球体,其自然表面是一个极其复杂且不规则的曲面。

地球形状更为精确的定义,尤其在通过同一大地水准面准定的高度时。 这一过程相当复杂的物理计算由地球质量产生。 地球上从一个区域到另一区域会有不同的厚度,因此会有不同的重力,而这重力影响了地球的形状。 所以,大地水准面表现了地球的重力场。 看一下上图,我们会发现地球的形状实际上是“扭曲”了的。 为了避免数学上的复杂性,在GIS通常将地球的形状显示为椭球。

大地测量中用水准测量方法得到的地面上各点的高程是依据一个理想的水准面来确定的, 该水准面通常称为“大地基准面”。 大地水准面是假定海水处于“完全”静止状态,将海水面延伸到大陆之下形成包围整个地球的连续表面; 大地水准面所包围的球体称为大地球体。 由于大地水准面上任何一点的铅垂线都与大地水准面成正交, 而铅垂线的方向又受地球内部质量分布不均匀的影响而有微小变化, 导致大地水准而产生微小的起伏。 因此,大地球体仍然是一个表面起伏的不规则球体,还不能直接作为投影的依据。

大地水准面是对地球形状的一个抽象、简化的表达。

地球椭球体

大地水准面所包围的形体,叫大地球体。大地水准面形状十分复杂,但从整体来看,起伏是微小的,它是一个很接近于绕自转轴(短轴) 旋转的椭球体。为了描述和表达地球表面,必须选择一个与地球形状、大小相接近的椭球体来近似代替它。 所以在测量和制图中就用旋转椭球来代替大地球体,这个旋转球体通常称地球椭球体,称椭球体。

将地球简化为球体,不适合创建比例尺大于1:2百万的地图。 旋转椭球和其它椭球试着调整地球的复杂形状,尽可以准确地满足数学计算。 因此,两极到地心的距离小于赤道到地心的距离。

地球椭球的两个主要参数为长轴半径a、短轴半径b,以及三个派生参数:

地球椭球体是对为了方便数学、计算机处理,而对大地水准面的进一步抽象、简化, 是一个数学表达。正是有了地球椭球体, 才使得目前阶段计算机处理地理空间数据有了可能。 因此,我们在以后所探讨的,只有地球椭球体这个概念。

当然也不排除一种可能性:未来有一天,测量技术与计算机技术发展到更高级的层面, 人们可以根据测量成果,在计算机内直接对大地水准面,甚至地球形状进行直接建模。

常用的地球椭球体参数

一个多世纪以来,各国学者与工程人员对地球进行了众多研究与测量,并提出了多组地球椭球体参数, 由于不同的地方,变形规律的不同, 因此各个国家根据本国的具体位置采用适合本国的投影的椭球体。

椭球模型有一个适用范围,它能够给地球上不同区域以最优的结果。 总之,我们是能够得到一个用于局部定位的、足够准确的基础。

我国使用的地球椭球体参数列表:

  1. 海福特椭球(1910) 我国52年以前基准椭球

  • a = 6378388m

  • b = 6356911.9461279m

  • alpha = 0.33670033670

  1. 克拉索夫斯基椭球(1940 Krassovsky)  北京54坐标系基准椭球

  • a = 6378245m

  • b = 6356863.018773m

  • alpha = 0.33523298692

  1. 1975年I.U.G.G推荐椭球(国际大地测量协会1975)  西安80坐标系基准椭球

  • a = 6378140m

  • b = 6356755.2881575m

  • alpha = 0.0033528131778

  1. WGS-84椭球(GPS全球定位系统椭球、17届国际大地测量协会) WGS-84 GPS 基准椭球

  • a = 6378137m

  • b = 6356752.3142451m

  • alpha = 0.00335281006247

注意

椭球体与基准面这两个概念是有区别的。 地图学上对地球上的抽象,第一次抽象为水准面(等重力面), 第二次抽象为椭球体(ellipsoid), 第三次抽象现在我认为是将椭球体进行定位之后, 所确定的具有明确的方向的椭球体, 它的要求能够很好的为当地区的地图制作服务。

在同一基准面上,基于数学公式的代数计算结果进行投影转换是可行的。 但是在不同基准面之间,是无法使用统一的数学公式计算来进行投影转换的。 基准面上空间位置的确定是测绘工作的成果, 不同基准面之间空间位置的对应关系转换必须由给定的点进行求算才可以。