1

I have an array of gridded data in XYZ format. I am trying to write shapefile (*.shp) file with multiple polygons with the pyshp Python libraries.

I am able to make it for only one polygon with the following code:

import shapefile


w = shapefile.Writer(shapeType=shapefile.POLYGONZ)
w.poly([[ [398010.0 7541990.0 280.4], [398010.0 7541980.0 281.5], [398020.0 7541980.0 280.9], [398020.0 7541990.0 279.8], [398010.0 7541990.0 280.4] ]], shapeType=15 )
w.field('NAME')
w.record('PolyZTest')
w.save('MyPolyZ')

However, I am not able for several polygons.

How can be done? For example, for the following four polygons?

POLYGONZ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))
POLYGONZ((398010.0 7542000.0 279.4, 398010.0 7541990.0 280.4, 398020.0 7541990.0 279.8, 398020.0 7542000.0 278.8, 398010.0 7542000.0 279.4))
POLYGONZ((398000.0 7541990.0 281.0, 398000.0 7541980.0 282.1, 398010.0 7541980.0 281.5, 398010.0 7541990.0 280.4, 398000.0 7541990.0 281.0))
POLYGONZ((398010.0 7541990.0 280.4, 398010.0 7541980.0 281.5, 398020.0 7541980.0 280.9, 398020.0 7541990.0 279.8, 398010.0 7541990.0 280.4))

If there are other easier command line options writing all polygons to a file, also could be a good choice.

Maybe the only possibility is to create shapefiles and then merge them?

  • have you looked up mpolygon – Brad Nesom Mar 3 '15 at 14:51
  • Brad, that is only for AutoCAD as I can see, deosn't it? I cannot find anything about mpolygon related with pyshp or similar. Can you explain what do you refer to? – iblasi Mar 3 '15 at 15:20
  • Do you refer to PolygonM? That rings do not allow XYZ space, only XYZ. – iblasi Mar 3 '15 at 17:56
3

If you want to use the WKT format, the correct syntax is:

"POLYGON Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))"

and not POLYGONZ( but it is very easy to transform your original format to correct WKT

poly = 'POLYGONZ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))'
poly.replace('POLYGONZ','POLYGON Z ')
print poly
'POLYGON Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))'

Then you can use any Python module which allows to convert the WKT format to valid coordinates for PyShp (Shapely, PyGeoif, ...)

With PyGeoif (easier to use):

from pygeoif import geometry
poly = 'POLYGON Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))' 
geom = geometry.from_wkt(poly)
print poly.exterior.coords
((398000.0, 7542000.0, 279.89999999999998), (398000.0, 7541990.0, 281.0), (398010.0, 7541990.0, 280.39999999999998), (398010.0, 7542000.0, 279.39999999999998), (398000.0, 7542000.0, 279.89999999999998))

So, first build a list with all the geometries

geometries = ['POLYGON Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))','POLYGONZ((398010.0 7542000.0 279.4, 398010.0 7541990.0 280.4, 398020.0 7541990.0 279.8, 398020.0 7542000.0 278.8, 398010.0 7542000.0 279.4))','POLYGONZ((398000.0 7541990.0 281.0, 398000.0 7541980.0 282.1, 398010.0 7541980.0 281.5, 398010.0 7541990.0 280.4, 398000.0 7541990.0 281.0))','POLYGONZ((398010.0 7541990.0 280.4, 398010.0 7541980.0 281.5, 398020.0 7541980.0 280.9, 398020.0 7541990.0 279.8, 398010.0 7541990.0 280.4))']

Then create the shapefile with a simple for loop:

import shapefile
w = shapefile.Writer(shapefile.POLYGONZ)
w.field('test','N')
for i,poly in enumerate(geometries):
      w.poly([geometry.from_wkt(poly).exterior.coords])
      w.record(i)

w.save('testpysh')

Result:

enter image description here

But is is as easy with Fiona and the geo_interface protocol ( GeoJSON-like)

import fiona
from pygeoif import geometry
# schema of the shapefile
schema = {'geometry': '3D Polygon','properties': {'test': 'int'}}
# creation of the shapefile
with fiona.open('testpysh2.shp','w','ESRI Shapefile', schema) as output:
    for i,poly in enumerate(geometries):
       output.write({'geometry':geometry.from_wkt(poly).__geo_interface__, 'properties':{'test':i}})

Control

for shape in fiona.open('testpysh2.shp):
    print shape
{'geometry': {'type': 'Polygon', 'coordinates': [[(398000.0, 7542000.0, 279.89999999999998), (398010.0, 7542000.0, 279.39999999999998), (398010.0, 7541990.0, 280.39999999999998), (398000.0, 7541990.0, 281.0), (398000.0, 7542000.0, 279.89999999999998)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'test', 0)])}
{'geometry': {'type': 'Polygon', 'coordinates': [[(398010.0, 7542000.0, 279.39999999999998), (398020.0, 7542000.0, 278.80000000000001), (398020.0, 7541990.0, 279.80000000000001), (398010.0, 7541990.0, 280.39999999999998), (398010.0, 7542000.0, 279.39999999999998)]]}, 'type': 'Feature', 'id': '1', 'properties': OrderedDict([(u'test', 1)])}
{'geometry': {'type': 'Polygon', 'coordinates': [[(398000.0, 7541990.0, 281.0), (398010.0, 7541990.0, 280.39999999999998), (398010.0, 7541980.0, 281.5), (398000.0, 7541980.0, 282.10000000000002), (398000.0, 7541990.0, 281.0)]]}, 'type': 'Feature', 'id': '2', 'properties': OrderedDict([(u'test', 2)])}
{'geometry': {'type': 'Polygon', 'coordinates': [[(398010.0, 7541990.0, 280.39999999999998), (398020.0, 7541990.0, 279.80000000000001), (398020.0, 7541980.0, 280.89999999999998), (398010.0, 7541980.0, 281.5), (398010.0, 7541990.0, 280.39999999999998)]]}, 'type': 'Feature', 'id': '3', 'properties': OrderedDict([(u'test', 3)])}
3

You are doing good so far. Just try to change the syntax of your polygonz declaration. I did it for you. Here is the complete working code with comments where neccessary:

# import shapefile library  
import shapefile

#create an in memory polygon shapefile
w = shapefile.Writer(shapeType=shapefile.POLYGONZ)
#add a name field to it
w.field('NAME')

#fill name field with polygon1 value
w.record('Polygon1')
#fill geometry of the feature with this polygon coordinates
w.poly([[[398000.0, 7542000.0, 279.9], [398000.0 ,7541990.0, 281.0], [398010.0 ,7541990.0, 280.4], [398010.0, 7542000.0, 279.4], [398000.0, 7542000.0, 279.9]]],shapeType=15)

w.record('Polygon2')
w.poly([[[398010.0, 7542000.0, 279.4], [398010.0 ,7541990.0, 280.4], [398020.0 ,7541990.0, 279.8], [398020.0, 7542000.0, 278.8], [398010.0, 7542000.0, 279.4]]],shapeType=15)

w.record('Polygon3')
w.poly([[[398000.0, 7541990.0, 281.0], [398000.0 ,7541980.0, 282.1], [398010.0 ,7541980.0, 281.5], [398010.0, 7541990.0, 280.4], [398000.0, 7541990.0, 281.0]]],shapeType=15)

w.record('Polygon4')
w.poly([[[398010.0, 7541990.0, 280.4], [398010.0 ,7541980.0, 281.5], [398020.0 ,7541980.0, 280.9], [398020.0, 7541990.0, 279.8], [398010.0, 7541990.0, 280.4]]],shapeType=15)

#save the save file to a file
w.save(r'D:\Cache\MyPolyZ')

I hope it helps

  • This answer is totally correct, but I found that the answer of @gene is more general in case of having WKT files or a file with several coordinates. Thank you for the answer! – iblasi Mar 5 '15 at 10:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.