python gis 经纬度 库_入门-Python-GIS坐标转换

前言

做GIS数据处理的同仁，不可避免的都会遇到坐标转换的问题，也许很多人遇到该问题，马上会使用各类GIS坐标转换的工具软件，甚至是GIS平台，比如ArcGIS，其实除非代转数据是未知坐标系(必须通过控制点进行配准)，只要是已知坐标系，都可采用proj4的开源实现来完成批处理转换，本文即以Python+pyproj来阐述如何进行批处理坐标转换。

Python

开始之前，提几句Python这门语言。这门语言像网红，一夜之间突然变成了全民语言，除了本职程序员，在全民学编程背景下首选都会选择Python，为什么？Python再次变成一等公民的原因，很大的功劳要算在近几年日益成熟的机器学习和深度学习平台，很多平台(比如TensorFlow)都会支持Python来做开发，而机器学习和人工智能又是未来热点，因此越来越多非科班的编程人员投身到这个行业，自然而然就会导致Python的兴起；

Python是非常容易入门和学习的，以TensorFlow为例，当它支持C++和Python的情况下，二选一，对于一个非科班人士，答案显而易见；

Python的出身是什么，就是主打数据批处理，而这点是提高生产效率的关键，除了IT，任何传统行业在日常工作中越来越多的需要进行数据处理尤其是批处理，当Excel完成不了时，也许Python是可能的出路。注意：Python并非唯一出路，例如坐标转换，当然也有java，.net，javascript(proj4.js)对应的proj4类库实现，这里不多赘述。

环境搭建

如果是从无到有的搭建步骤如下：Python安装，推荐Python3以上，当前Python37(与Python27，有大量不兼容的函数和API，注意ArcGIS10.X平台安装，默认会安装Python27，但不会冲突)官网

IDE环境安装，推荐VS Code(Free)，Pycharm(Buy)。

在VS Code下，安装Python扩展。

设置VS Code对应的Python编译器，此步在安装有多个Python版本编译器时，必须。

安装pyproj类库，pip3 install pyproj。

应用场景

为了说明如何利用proj4来完成批处理转换，暂将场景设置如下：

记录一组当地坐标系的坐标的文本文件(此处暂考虑文本文件，其实只要是有格式说明的或白皮书的GIS格式，都可以采用批处理来完成，只不过添加相应的格式读取类库来进行数据预处理，比如shp，geojson等等，选择文本文件的原因，是本文关注点是坐标转换。)，如何将这组坐标叠加到高德地图上？(高德地图其实是web mercator，但按国测局要求进行了偏移，网络上大家称为国测局gcj02)

坐标转换流程

地方坐标系->WGS84->WGS84偏移(GCJ02 经纬度)->Web Mercator偏移(GCJ02 投影后平面坐标，这步其实可选)

脚本代码

import os

from pyproj import CRS

from pyproj import Transformer

from converter import wgs84_to_gcj02 #参见注意事项

input_file = './input.txt'

output_file = './output.txt'

#当地坐标系转WGS84

from_crs = CRS.from_wkt('PROJCS["local",GEOGCS["GCS_Xian_1980",DATUM["D_Xian_1980",SPHEROID["Xian_1980",6378140.0,298.257]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",135005.0014],PARAMETER["False_Northing",-1999781.9795],PARAMETER["Central_Meridian",109.75],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]')

to_crs = CRS.from_epsg(4326)

transformer = Transformer.from_crs(from_crs, to_crs, always_xy=True)

#WGS84转Web Mercator

from_crs_2 = CRS.from_epsg(4326)

to_crs_2 = CRS.from_epsg(3857)

transformer_2 = Transformer.from_crs(from_crs_2, to_crs_2, always_xy=True)

with open(output_file, "w") as fo:

with open(input_file, "r") as fi:

while True:

line = fi.readline() # 逐行读取

if not line:

break

else:

array = line.split(",") # x,y 逗号分隔

x1,y1 = transformer.transform(array[0], array[1]) # 当地坐标系转WGS84

x2,y2 = wgs84_to_gcj02(x1, y1) # gcj02 坐标偏移

x3,y3 = transformer_2.transform(x2, y2) # WGS84转Web Mercator

fo.write(",".join(["{:.6f}".format(x3),"{:.6f}".format(y3),'\n'])) # 输出到新文件

print('All Done!')

注意：此处借用github上WGS84转GCJ02的Python脚本，请自行下载。

定义坐标系

常用两种方式：

API: CRS.from_wkt

wkt坐标系描述字符串，适合自定义，也可来源于prj文件。

API: CRS.from_epsg

epsg请查询epsg官网，可以认为官方给通用坐标系颁发的一个唯一编码，请记住WGS84为4326，Web Mercator为3857。