Shp文件坐标转化

背景

当你的路网数据和瓦片地图有偏差时候，怎么办呢?我们肯定能想到：把坐标转化到同一个坐标系里面就好啦。没错，原理是这样的。

但是，我们应该怎么把Shp文件里面的坐标点提取出来，然后转化呢？
这就需要用到Python的shapefile库了。

示例文件：链接：https://pan.baidu.com/s/1hbxomqLjdy3X_Kx4WsnHxQ
提取码：8z5a

PyShp官方文档：https://pypi.org/project/pyshp/

下面，我用一个有偏差的Shp文件，把文件的经纬度全部提取出来，然后经过转化后重新生成一个Shp文件。

代码实现

主代码：

import shapefile
# import untitled.ConvertLag_Lat as CVLL
import ConvertLag_Lat as CVLL

# 读取文件的路径
readPath = r"示例文件路径\城市快速路.shp"
# 要写入新的文件存放路径
writeFilePath = r"你想输出的路径\纠偏后_城市快速路.shp"




'''读取文件的操作'''
def ShpReadFile(filePath):
    sf = shapefile.Reader(filePath,encoding="gbk")
    # sf = shapefile.Reader(filePath)
    # 以ShapeRecords格式读入(既可以访问shape几何对象，又可以访问record属性)
    ShapeRecords = sf.shapeRecords()
    return ShapeRecords,sf



'''从读取的文件把相关的内容写入到新的文件中'''
if __name__ == '__main__':
    # 得到读取文件的对象
    ShapeRecords,sf = ShpReadFile(readPath)
    Shapes = sf.shapes()

    # 写入文件的操作,并且是以线的方式
    w = shapefile.Writer(writeFilePath)
    # 属性列直接赋值
    w.fields = sf.fields[1:]
    # print(sf.fields[1:])
    # **自动平衡功能
    w.autoBalance = 1
    # 遍历所有的记录
    for shaperec in sf.iterShapeRecords():
        # 属性赋值
        w.record(*shaperec.record)
        # 先取出每个几何对象的点的集合
        sub_points = shaperec.shape.points
        # 使用geo格式查看
        # print("parts："+str(shaperec.__geo_interface__))

        # 点的总集合
        poly_list = []
        # 遍历点的集合
        for point in sub_points:
            # 把WGS84经纬度转化为GCJ02
            rig_Point = CVLL.wgs84togcj02(point[0],point[1])
            poly_coord = [float(rig_Point[0]),float(rig_Point[1])]
            poly_list.append(poly_coord)
        # 把转化后的点添加进去(画线)
        w.line([poly_list])
    # 关闭写操作
    w.close()

经纬度代码：

import math
x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626  # π
a = 6378245.0  # 长半轴
ee = 0.00669342162296594323  # 扁率


# Web墨卡托转经纬度
def Mercator2lonLat(mercatorX,mercatorY):
    x = mercatorX/20037508.34*180
    y = mercatorY/20037508.34*180
    y= 180/math.pi*(2*math.atan(math.exp(y*math.pi/180))-math.pi/2)
    return [x,y]


# WGS84转GCJ02
def wgs84togcj02(lng, lat):
    """
    WGS84转GCJ02(火星坐标系)
    :param lng:WGS84坐标系的经度
    :param lat:WGS84坐标系的纬度
    :return:
    """
    if out_of_china(lng, lat):  # 判断是否在国内
        return lng, lat
    dlat = transformlat(lng - 105.0, lat - 35.0)
    dlng = transformlng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * pi
    magic = math.sin(radlat)
    magic = 1 - ee * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
    dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return [mglng, mglat]

def out_of_china(lng, lat):
    """
    判断是否在国内，不在国内不做偏移
    :param lng:
    :param lat:
    :return:
    """
    if lng < 72.004 or lng > 137.8347:
        return True
    if lat < 0.8293 or lat > 55.8271:
        return True
    return False


def transformlat(lng, lat):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
        0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * pi) + 40.0 *
            math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 *
            math.sin(lat * pi / 30.0)) * 2.0 / 3.0
    return ret


def transformlng(lng, lat):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
        0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
    ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 *
            math.sin(2.0 * lng * pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * pi) + 40.0 *
            math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 *
            math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return ret

效果图

Shp文件坐标转化

红色的为经过纠偏后的Shp文件，很明显，纠偏后就和瓦片地图重叠了。
再看看Shp文件数据，看看有没有损坏？

Shp文件坐标转化

OK，数据完好无损。大功告成！