Multiple polygon to shapefile in XYZ format

Question

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?

Answer 1

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:

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)])}

Answer 2

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