[GIS] Multiple polygon to shapefile in XYZ format

polygonpyshpshapefile

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?

Best Answer

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