pyhton 创建shp文件并投影，批量进行坐标转换与添加投影信息，合并shp文件，分割shp文件

代码如下，注释是经过百度翻译的中文。这是我的python地理分析课程作业之一。

如何获取投影信息，可以通过arcgis输出投影信息，更换投影信息的时候主要也要进行坐标转换，尤其从地理坐标系转到投影坐标系。

import shapefile
import utm


# Add projection information, the spatial coordinate system is wgs-1984-mercator projection CGS
def projection(filename, ind):
    with open(filename, 'w') as ww:
        if ind == 1:
            ww.write('GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],'
                     'PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]')
        if ind == 2:
            ww.write('PROJCS["WGS_1984_UTM_Zone_50N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",'
                     '6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],'
                     'PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER['
                     '"False_Northing",0.0],PARAMETER["Central_Meridian",117.0],PARAMETER["Scale_Factor",0.9996],'
                     'PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]')
    with open(filename.replace(".prj", ".cpg"), 'w') as ww1:
        ww1.write("UTF-8")


# Create a series of points using pyshp library. Only geometric figures of the same type can be merged.
# Shapetype=1 represents points
with shapefile.Writer("./data/school", shapeType=1, encoding='UTF-8') as w:
    # Define two attributes, namely the name and number of the point
    w.field('name', 'C', 50, 0)
    w.field('id', 'N', 10, 0)
    # Storing geometric locations and attributes
    w.point(117.239252, 31.510307)
    w.record('合肥工业大学', 152453)
    w.point(117.169232, 31.550335)
    w.record('安徽大学', 312854)
    w.point(117.288322, 31.821104)
    w.record('三联学院', 323664)
    projection("./data/school.prj", 1)

with shapefile.Writer("./data/hospital", shapeType=1, encoding='UTF-8') as w:
    w.field('name', 'C', 50, 0)
    w.field('id', 'N', 10, 0)
    w.point(117.262517, 31.841903)
    w.record('安徽省人民医院', 135978)
    w.point(117.260985, 31.861645)
    w.record('安徽省第二医院', 132043)
    w.point(117.269712, 31.863475)
    w.record('工大附属医院', 343124)
    projection("./data/hospital.prj", 1)

with shapefile.Writer("./data/market", shapeType=1, encoding='UTF-8') as w:
    w.field('name', 'C', 50, 0)
    w.field('id', 'N', 10, 0)
    w.point(117.244368, 31.846191)
    w.record('万达购物广场', 136852)
    w.point(117.297992, 31.836775)
    w.record('家乐购', 246951)
    w.point(117.206234, 31.870627)
    w.record('直达便利店', 133456)
    projection("./data/market.prj", 1)

# Perform batch projection
files = ["./data/hospital", "./data/school", "./data/market"]
for file in files:
    r = shapefile.Reader(file)
    w = shapefile.Writer(file + 'UTM-50N', encoding='UTF-8')
    w.fields = r.fields
    for i in r.iterShapeRecords():
        lon, lat = i.shape.points[0]
        y, x, zone, band = utm.from_latlon(lat, lon)
        w.point(x, y)
        w.record(*i.record)
    w.close()
    projection(file + 'UTM-50N' + '.prj', 2)

# Merge three point layers into city layer
with shapefile.Writer("./city", shapeType=1, encoding='UTF-8') as w:
    for m in files:
        r = shapefile.Reader(m + 'UTM-50N')
        w.fields = r.fields
        for i in r.iterShapeRecords():
            w.shape(i.shape)
            w.record(*i.record)
    projection('city.prj', 2)

# Divide the city layer and save points with X coordinate greater than 4080694
with shapefile.Writer('./city_spilt', encoding='UTF-8') as w:
    r = shapefile.Reader('city')
    w.fields = r.fields
    for i in r.iterShapeRecords():
        if i.shape.points[0][0] > 4080694:
            w.shape(i.shape)
            w.record(*i.record)
    projection('./city_spilt.prj', 2)