读面状shp的中心经纬度，并增加shp字段信息，把经纬度写进shp属性表里

#!/usr/bin/python
# -*- coding: utf-8 -*-

import shapefile
from osgeo import gdal, ogr, osr
import numpy as np
from sympy import *
import math
from math import cos
from mpl_toolkits.basemap import Basemap
import matplotlib
import matplotlib.pyplot as plt
import os
import sys

file = r'1.shp'
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(file, 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("LON", ogr.OFTReal)#####在shp里增加一列经度坐标
new_field.SetWidth(32)####设置长度
new_field.SetPrecision(7)####设置小数点位数
layer.CreateField(new_field)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("LAT", ogr.OFTReal)
new_field.SetWidth(32)
new_field.SetPrecision(7)  # 设置小数点
layer.CreateField(new_field)
file=shapefile.Reader(r'1.shp')
shapes=file.shapes()
records=file.records()
# print(records)


for i in range(len(shapes)):
    points = shapes[i].points  #各个点的经纬度
    lon = []
    # print(lon)
    lat = []
    # a = np.array(points, dtype=np.float32)#列表转二维数组
    # a1 = a[ : , 0 ]#提取二维数组中的第一列
    # print(a1)
    [lon.append(points[i][0]) for i in range(len(points))]#表示在原来的列表末尾追加新的对象,并且提取出lon的列表
    [lat.append(points[i][1]) for i in range(len(points))]
    # print(lon)

    lon = np.array(lon).reshape(-1, 1)
    # print(lon)
    lon1 = np.max(lon)
    lon2 = np.min(lon)
    lon3 = (lon1 +lon2)/2
    # print(lon3)

    lat = np.array(lat).reshape(-1, 1)
    # print(lat)
    lat1 = np.max(lat)
    lat2 = np.min(lat)
    lat3 = (lat1 +lat2)/2
    # print(lat3)
    feature = layer[i]
    feature.SetField("LON", lon3)  # 将经纬度添加到属性表中
    feature.SetField("LAT", lat3)
    layer.SetFeature(feature)
dataSource = None

作者在线提问，有没有人用过w=shapefile.Writer(r’1.shp’)这个函数，本来想用这个函数写经纬度，以及shp的，结果发现投影坐标丢失。。。然后给他加上投影坐标或地理坐标后。shp的地理位置就跑走不对了，，，，请问有没有人知道怎么破