pyshp写投影信息

pyshp库对shp文件的操作本身不包含投影信息，如果需要增加投影信息的话需要单独写.prj文件

1、使用pyshp写shp文件（包含shp, dbf, shx）

"""
using pyshp create shpfile
(1) type:You can reference shape types by the numbers or by constants defined by PyShp:
shapefile.NULL = 0 shapefile.POINT = 1 shapefile.POLYLINE = 3 shapefile.POLYGON = 5
shapefile.MULTIPOINT = 8 shapefile.POINTZ = 11 shapefile.POLYLINEZ = 13 shapefile.POLYGONZ = 15
shapefile.MULTIPOINTZ = 18 shapefile.POINTM = 21 shapefile.POLYLINEM = 23 shapefile.POLYGONM = 25
shapefile.MULTIPOINTM = 28 shapefile.MULTIPATCH = 31'''
(2) field setting
Field name: the name describing the data at this column index.
Field type: the type of data at this column index. Types can be: Character, Numbers, Longs, Dates, or Memo. The “Memo” type has no meaning within a GIS and is part of the xbase spec instead.
Field length: the length of the data found at this column index. Older GIS software may truncate this length to 8 or 11 characters for “Character” fields.
Decimal length: the number of decimal places found in “Number” fields.

"""

import shapefile

outshp = r'D:\za\testshp\a.shp'

w = shapefile.Writer(shapeType=5)

#设置字段，最大长度为254，C为字符串
w.field('FIRST_FLD')
w.field('SECOND_FLD','C','40')
#添加几何和添加字段信息，添加两个示例，字段顺序区分先后
w.poly([[[122,37], [117,36], [115,32],[118,20], [123,24],[122,37]]])
w.record('First','Point')
w.poly([[[123,37], [118,36], [116,32],[119,20], [124,24],[123,37]]])
w.record('Second','Point')
#保存
w.save(outshp)

2、写投影信息（.prj文件）

通过.prj文件设置投影，需要写入一个wkt字符串。什么是WKT呢？英文是well-known text，即用“大家都知道的文本”描述投影信息，本质是一串字符串。这英文着实有点大白话… 例如

wkt="""PROJCS["WGS_1984_UTM_zone_50N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32650"]]'"""

1）来自于定义的投影wkt字符串

import osr

##gdal的GetProjection()返回的是wkt字符串，需要ImportFromWkt
proj = osr.SpatialReference()
proj.ImportFromEPSG(4326)
#或 proj.ImportFromProj4(proj4str)等其他的来源
wkt = proj.ExportToWkt()

#写出prj文件
f = open(outshp.replace(".shp",".prj"), 'w')
f.write(wkt)
f.close()

2）来自于已有shp文件

in_shp_file = 'D:/shape_of_input.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open(in_shp_file)
# from Layer
layer = dataset.GetLayer()
spatialRef = layer.GetSpatialRef()
wkt = spatialRef.ExportToWkt()

#写出prj文件
f = open(outshp.replace(".shp",".prj"), 'w')
f.write(wkt)
f.close()

3）来自于影像

import gdal

image_path = 'D:/input_image_path.tif'
raster = gdal.Open(image_path)
wkt = raster.GetProjection()

#写出prj文件
f = open(outshp.replace(".shp",".prj"), 'w')
f.write(wkt)
f.close()

参考：

[1] pyshp创建面shp文件并设置投影

[2] map-projections

[3] pyshp - CreatePRJfiles.wiki

Creating a .prj file

The shapefile format does not allow for specifying the map projection of the data. As a workaround ESRI invented another file extension called ".prj" which is short for projection. The prj file contains a WKT (Well-Known Text) string which has all the parameters for the map projection. So the format is quite simple and was created by the Open GIS Consortium.

But there are several thousand "commonly-used" map projections which were standardized by the European Survey Petroleum Group (EPSG). And there's no way to accurately detect the projection from the coordinates in the shapefile. For these reasons the Python Shapefile Library does not currently handle prj files.

If you need a prj file the easiest thing to do is write one yourself. The following example creates a simple point shapefile and then the corresponding prj file using the WGS 84 "unprojected" WKT.

I've thought about adding the ability to optionally write prj files but the list of "commonly-used" WKT strings is over .5 megs and would be bigger than the shapefile library itself.

The easiest thing to do right now is just figure out what WKT string you need for your data and quickly write a file:

``` import shapefile as sf filename = 'test/point'

create the shapefile

w = sf.Writer(sf.POINT) w.point(37.7793, -122.4192) w.field('FIRST_FLD') w.record('First','Point') w.save(filename)

create the PRJ file

prj = open("%s.prj" % filename, "w") epsg = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]' prj.write(epsg) prj.close() ```