1、地理坐標(biāo)系與投影坐標(biāo)系

空間參考中主要包含大地水準(zhǔn)面、地球橢球體、投影坐標(biāo)系等幾部分內(nèi)容。地圖投影就是把地球表面的任意點(diǎn)，利用一定數(shù)學(xué)法則，轉(zhuǎn)換到地圖平面上的理論和方法，一般有兩種坐標(biāo)系來進(jìn)行表示，分別是地理坐標(biāo)系和投影坐標(biāo)系。如下圖所示，描述了地理坐標(biāo)系與投影坐標(biāo)系之間的關(guān)系。

圖1 地理坐標(biāo)與投影坐標(biāo)關(guān)系

地理坐標(biāo)系即球體坐標(biāo)，投影坐標(biāo)即投影后的平面坐標(biāo)；常見的地理坐標(biāo)GPS坐標(biāo)，常見的投影坐標(biāo)墨卡托投影。

常用的坐標(biāo)系為地理坐標(biāo)系（Geograpic Coordinate System,簡稱GCS）和投影坐標(biāo)系（Projected Coordinate System，簡稱PCS）。

1.1、地理坐標(biāo)系統(tǒng)

地理坐標(biāo)系統(tǒng)（GCS）用一個三維的球面來確定地物在地球上的位置，地面點(diǎn)的地理坐標(biāo)有經(jīng)度、緯度、高程構(gòu)成。地理坐標(biāo)系統(tǒng)與選擇的地球橢球體和大地基準(zhǔn)面有關(guān)。橢球體定義了地球的形狀，而大地基準(zhǔn)面確定了橢球體的中心。

　　地理坐標(biāo)系 (GCS) 使用三維球面來定義地球上的位置。GCS中的重要參數(shù)包括角度測量單位、本初子午線和基準(zhǔn)面（基于旋轉(zhuǎn)橢球體）。地理坐標(biāo)系統(tǒng)中用經(jīng)緯度來確定球面上的點(diǎn)位，經(jīng)度和緯度是從地心到地球表面上某點(diǎn)的測量角。球面系統(tǒng)中的水平線是等緯度線或緯線，垂直線是等經(jīng)度線或經(jīng)線。這些線包絡(luò)著地球，構(gòu)成了一個稱為經(jīng)緯網(wǎng)的格網(wǎng)化網(wǎng)絡(luò)。

GCS中經(jīng)度和緯度值以十進(jìn)制度為單位或以度、分和秒 (DMS) 為單位進(jìn)行測量。緯度值相對于赤道進(jìn)行測量，其范圍是 -90°（南極點(diǎn)）到 +90°（北極點(diǎn)）。經(jīng)度值相對于本初子午線進(jìn)行測量。其范圍是 -180°（向西行進(jìn)時(shí)）到 180°（向東行進(jìn)時(shí)）。

1.2、投影坐標(biāo)系統(tǒng)

投影坐標(biāo)系統(tǒng)是根據(jù)某種映射關(guān)系，將地理坐標(biāo)系統(tǒng)中由經(jīng)緯度確定的三維球面坐標(biāo)投影到二維的平面上所使用的坐標(biāo)系統(tǒng)。在該坐標(biāo)系統(tǒng)中，點(diǎn)的位置是由(x,y,z)坐標(biāo)來確定的。由于投影坐標(biāo)是將球面展會在平面上，因此不可避免會產(chǎn)生變形。這些變形包括3種：長度變形、角度變形以及面積變形。通常情況下投影轉(zhuǎn)換都是在保證某種特性不變的情況下犧牲其他屬性。根據(jù)變形的性質(zhì)可分為等角投影、等面積投影等。

我國的基本比例尺地形圖（1：5千，1：1萬，1：2.5萬，1：10萬，1：25萬，1：50萬，1：100萬）中，大于或等于1：50萬均采用高斯-克呂格投影（Gauss_Kruger），又叫橫軸墨卡托投影（Transverse Mercator）；1：100萬的地形圖采用正軸等角圓錐投影，又叫蘭勃特投影（Lambert Conformal Conic）；海上小于50萬的地形圖多用正軸等角圓柱投影，又叫墨卡托投影（Mercator）。在開發(fā)GIS系統(tǒng)中應(yīng)該采用與我國基本比例尺地形圖系列一致的地圖投影系統(tǒng)。

2、坐標(biāo)系轉(zhuǎn)換

根據(jù)參考文章，封裝成一個類使用。

import numpy as np
from osgeo import gdal, osrclass GdalTif(object):def __init__(self, tif_path):self.tif_path = tif_pathself.dataset = self._read_tif(tif_path)@staticmethoddef _read_tif(tif_path):""" 讀取GDAL地理數(shù)據(jù):param tif_path: tif圖像路徑:return: GDAL地理數(shù)據(jù)"""dataset = gdal.Open(tif_path)if dataset is None:print(tif_path + " 文件無法打開")return datasetdef get_geo_trans(self):""" 獲取仿射矩陣信息:return: 仿射矩陣信息有六個參數(shù)，描述的是柵格行列號和地理坐標(biāo)之間的關(guān)系:(0：左上角橫坐標(biāo)（投影坐標(biāo)，經(jīng)度）；1：像元寬度,影像東西/水平方向分辨率；2：行旋轉(zhuǎn)，如果圖像北方朝上，該值為0；3：左上角縱坐標(biāo)（投影坐標(biāo)，緯度）；4：列旋轉(zhuǎn)，如果圖像北方朝上，該值為0；5：像元高度,影像南北/垂直方向分辨率；)如果圖像不含地理坐標(biāo)信息，默認(rèn)返回值是：(0,1,0,0,0,1)"""return self.dataset.GetGeoTransform()def get_projection(self):""" 獲取數(shù)據(jù)投影信息:return: INFO"""return self.dataset.GetProjection()def get_srs_pair(self):""" 獲得給定數(shù)據(jù)的投影參考系和地理參考系:param dataset: GDAL地理數(shù)據(jù):return: 投影參考系和地理參考系"""prosrs = osr.SpatialReference()prosrs.ImportFromWkt(self.dataset.GetProjection())geosrs = prosrs.CloneGeogCS()return prosrs, geosrsdef imagexy2geo(self, row, col):""" 根據(jù)GDAL的六參數(shù)模型將影像圖上坐標(biāo)（行列號）轉(zhuǎn)為投影坐標(biāo)或地理坐標(biāo)（根據(jù)具體數(shù)據(jù)的坐標(biāo)系統(tǒng)轉(zhuǎn)換）:param dataset: GDAL地理數(shù)據(jù):param row: 像素的行號(i):param col: 像素的列號(j):return: 行列號(row, col)對應(yīng)的投影坐標(biāo)或地理坐標(biāo)(x, y) - WGS84"""trans = self.dataset.GetGeoTransform()x = trans[0] + col * trans[1] + row * trans[2]y = trans[3] + col * trans[4] + row * trans[5]return (x, y)def geo2imagexy(self, x, y):""" 根據(jù)GDAL的六 參數(shù)模型將給定的投影或地理坐標(biāo)轉(zhuǎn)為影像圖上坐標(biāo)（行列號）:param x: 投影或地理坐標(biāo)x:param y: 投影或地理坐標(biāo)y:return: 影坐標(biāo)或地理坐標(biāo)(x, y)對應(yīng)的影像圖上行列號(row, col)"""trans = self.dataset.GetGeoTransform()a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])b = np.array([x - trans[0], y - trans[3]])return np.linalg.solve(a, b)[::-1]  # 使用numpy的linalg.solve進(jìn)行二元一次方程的求解def geo_to_lonlat(self, x, y):""" 將投影坐標(biāo)轉(zhuǎn)為經(jīng)緯度坐標(biāo)（具體的投影坐標(biāo)系由給定數(shù)據(jù)確定）:param x: 投影坐標(biāo)x:param y: 投影坐標(biāo)y:return: 投影坐標(biāo)(x, y)對應(yīng)的經(jīng)緯度坐標(biāo)(lon, lat)"""prosrs, geosrs = self.get_srs_pair()ct = osr.CoordinateTransformation(prosrs, geosrs)coords = ct.TransformPoint(x, y)return coords[:2]def lonlat2geo(self, lon, lat):""" 將經(jīng)緯度坐標(biāo)轉(zhuǎn)為投影坐標(biāo)（具體的投影坐標(biāo)系由給定數(shù)據(jù)確定）:param dataset: GDAL地理數(shù)據(jù):param lon: 地理坐標(biāo)lon經(jīng)度:param lat: 地理坐標(biāo)lat緯度:return: 經(jīng)緯度坐標(biāo)(lon, lat)對應(yīng)的投影坐標(biāo)"""prosrs, geosrs = self.get_srs_pair()ct = osr.CoordinateTransformation(geosrs, prosrs)coords = ct.TransformPoint(lon, lat)return coords[:2]def main():tif_file = r"C:\Users\i\Desktop\dsm.tif"gda = GdalTif(tif_file)print("數(shù)據(jù)投影信息:", gda.get_projection())pt = (300, 500)print("pt = ", pt)print("圖上坐標(biāo) -> 投影坐標(biāo)：")point = gda.imagexy2geo(*pt)print("WGS84 point:", point)print("投影坐標(biāo) -> 圖上坐標(biāo)：")pt = gda.geo2imagexy(*point)print("pixel pt:", pt)print("投影坐標(biāo) -> 經(jīng)緯度：")coord = gda.geo_to_lonlat(*point)print("gps coord:", coord)print("經(jīng)緯度 -> 投影坐標(biāo)：")point = gda.lonlat2geo(*coord)print("point:", point)if __name__ == '__main__':main()

測試結(jié)果：

數(shù)據(jù)投影信息: PROJCS["CGCS2000 / 3-degree Gauss-Kruger CM 111E",GEOGCS["China Geodetic Coordinate System 2000",DATUM["China_2000",SPHEROID[
"CGCS2000",6378137,298.257222101,AUTHORITY["EPSG","1024"]],AUTHORITY["EPSG","1043"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["d
egree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4490"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_orig
in",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",1],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["
metre",1,AUTHORITY["EPSG","9001"]],AXIS["Northing",NORTH],AXIS["Easting",EAST],AUTHORITY["EPSG","4546"]]
pt =  (300, 500)
圖上坐標(biāo) -> 投影坐標(biāo)：
WGS84 point: (402657.7898284429, 2052764.045716705)
投影坐標(biāo) -> 圖上坐標(biāo)：
pixel pt: [300. 500.]
投影坐標(biāo) -> 經(jīng)緯度：
gps coord: (3.535317397929979, 124.8387146988415)
經(jīng)緯度 -> 投影坐標(biāo)：
point: (402657.78982844297, 2052764.0457167062)

參考文章：

1、談?wù)劦乩碜鴺?biāo)和投影坐標(biāo)

http://www.360doc.com/content/22/1214/01/65719465_1060178233.shtml